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INTRODUCTION 


This document describes progress in the development of finite element codes for the 
prediction of near and far field acoustic radiation from the inlet and aft fan ducts of turbofan 
engines. The report consists of nine papers which have appeared in archival journals and 
conference proceedings, or are presently in review for publication. Brief summaries of each 
paper are given here. 

1. Aft Fan Duct Acoustic Radiation. (Journal of Sound and Vibration 1998 213(2), 

235- 257) 

Details are given of a finite element code which has been developed for the prediction 
of the radiated acoustic field from the aft fan duct of a turbofan engine. A new technique 
based on a penalty method is introduced to enforce the condition of continuity of acoustic 
pressure across the shear layer which bounds the jet. 

2. Mapped Infinite Wave Envelope Elements for Acoustic Radiation in a Uniformly 

Moving Medium. (Journal of Sound and Vibration 1999 224(4), 665-687) 

Variable order mapped infinite wave envelope elements are introduced for finite 
element modeling of acoustic radiation in a uniformly moving medium. The elements are 
applied to the problem of turbofan inlet radiation, and are shown to provide an effective non- 
reflecting boundary condition which allows substantial reduction of the FEM mesh in the 
near field. Results are shown for the acoustic pressure in the near field. 

3. A Reflection Free Boundary Condition for Propagation in Uniform Flow Using Mapped 

Infinite Wave Envelope Elements. (Journal of Computational Acoustics 2000, 8(1), 

25-42.) 

Variable order mapped infinite wave envelope elements are introduced for finite 
element modeling of acoustic radiation in a uniformly moving medium. The elements are 
applied to the problem of turbofan inlet and aft fan duct radiation, and are shown to provide 
an effective non-reflecting boundary condition which allows substantial reduction of the 
FEM mesh in the near field. Results are shown for the acoustic pressure in the near field and 
in the far field. 



4. A Numerical Comparison Between Multiple-Scales and FEM Solution for Sound 

Propagation in Lined Flow Ducts. (AIAA Paper 99-1821, 1999 AIAA Aeroacoustics 

Conference, Seattle. Provisionally accepted for the Journal of Fluid Mechanics) 

An analytical solution based on the method of multiple-scales for acoustic 
propagation in nonuni form ducts with compressible potential mean flow is compared against 
a Finite element solution. This provides a useful benchmark for the FEM results in cases 
when scattering of the incident mode does not dominate the transmitted acoustic field. In 
general this will be true for incident modes which are nearly cut off. 

5. Acoustic Propagation at High Frequencies in Ducts. (AIAA Paper 2000-1953, 2000 

AIAA Aeroacoustics Conference, Maui, Hawaii) 

The problem of acoustic propagation in ducts at high non-dimensional frequencies 
is examined. It is found that in FEM models using quadratic elements, good solutions for the 
acoustic potential are achieved using the conventional 10 node per wave length rule of 
thumb. However, good solutions, via postprocessing, for acoustic pressure require 
substantially increased mesh density. Cubic and quartic elements are examined and it is 
found that cubic elements offer are more efficient than quadratic elements for acoustic 
potential and offer a substantial improvement in acoustic pressure post-processing. Non- 
dimensional frequencies up to 1 00 are considered. 

6. The Boundary Condition at an Impedance Wall in a Nonuniform Duct with Potential 

Flow. (Journal of Sound and Vibration, in review) 

The acoustic boundary condition at an impedance wall in a nonuniform duct with 
compressible mean flow is implemented in a weighted residuals finite element formulation. 
The boundary condition appears to require data which includes the tangential derivative of 
the tangential mean flow velocity, the normal derivative of the normal component of mean 
flow velocity, and the derivatives of the mean flow density and the boundary admittance 
along the boundary. It is shown that it can be substantially simplified to eliminate the 
tangential derivatives of mean flow properties and to eliminate completely the normal 
component of mean flow velocity. Implementation of the boundary condition is shown to 
involve no difficulty. 

7. A Reverse Flow Theorem and Acoustic Reciprocity in Compressible Potential Flows. 

(Journal of Sound and Vibration, in review) 

A reverse flow theorem for acoustic propagation in compressible potential flow has 
been obtained directly from the field equations without recourse to energy conservation 
arguments. A reciprocity theorem for the scattering matrix for propagation of acoustic modes 
in a duct with either acoustically rigid walls or acoustically absorbing walls follows. It is 
found that for a source at a specific end of the duct, suitably scaled reflection matrices in 
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direct and reverse flow have a reciprocal relationship. Scaled transmission matrices obtained 
for direct flow and reversed flow with simultaneous switching of source location from one 
end to the other also have a reciprocal relationship. 

Reciprocal relations provide an excellent benchmark for verification of acoustic propagation 
computations. Numerical verification of the reciprocal relationships is given in a companion 
paper. 

8. Reciprocity and Acoustic Power in One Dimensional Compressible Potential Flows. 

(Journal of Sound and Vibration, in review) 

A reverse flow theorem for one dimensional acoustic propagation in compressible 
potential flow has been obtained directly from the field equations without recourse to energy 
conservation arguments. Reciprocity relationships for the scattering coefficients for 
propagation are derived. It is found that for a source at a specific end of the duct, suitably 
scaled reflection coefficients in direct and reverse flow have a reciprocal relationship. Scaled 
transmission coefficients obtained for direct flow and reversed flow with simultaneous 
switching of source location from one end to the other also have a reciprocal relationship. 
Reciprocal relations and power conservation arguments are used to show that scaled power 
reflection and transmission coefficients are invariant to flow reversal and switching of source 
location from one end of the duct to the other. Numerical verification of the reciprocal 
relationships is given in a companion paper in which multiple mode propagation and one 
dimensional propagation are considered. 

9. Numerical Experiments on Acoustic Reciprocity in Compressible Potential Flows. 

(Journal of Sound and Vibration, in review) 

A reciprocity theorem for the scattering matrix for propagation of acoustic modes in 
a duct with acoustically hard walls or with acoustically absorbing walls has been given in a 
companion publication. It was found that for a source at a specified end of the duct, suitably 
scaled reflection matrices in direct and reverse flow have a reciprocal relationship. Scaled 
transmission matrices obtained for direct flow and reversed flow with simultaneous 
switching of source location from one end to the other also have a reciprocal relationship. 
A reverse flow theorem for the equivalent one dimensional propagation model, which is a 
good approximation to the three dimensional model at low frequencies, was also obtained. 
In this case, using reciprocity and acoustic power conservation arguments it is additionally 
found that the acoustic power transmission coefficient is the same for a source at either end 
of the duct for a given flow direction. This result leads to an invariance theorem which 
relates acoustic power propagated due to sources of equal pressure amplitude at the two ends 
of the duct. Numerical verification of these reciprocal relationships is given here for 
propagation in axially symmetric (circular and annular) ducts with multi-modal propagation 
and at low frequencies when a one dimensional model is appropriate. 
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AFT FAN DUCT ACOUSTIC RADIATION 


W. Eversman 


Department of Mechanical and Aerospace Engineering 
University of Missouri- Rolla, Rolla , MO 


and Engineering 
65401, U.S.A. 


Mechanics, 


AND 


D. Okunbor 

Department of Computer Science, University of Missouri- Rolla, Rolla . MO 63401, U.S.A 


( Received 13 June 1997, ami in final form 18 December 1997) 


A finite element code has been developed for the prediction of the radiated “tic field 
from the aft fan duct of a turbofan engine. The acoustic field is modelled based on the 
assumption that the steady flow in and around the nacelle is irrotational as is the acoustic 
perturbation. The geometry of the nacelle is ax, symmetric and the acoustic source ,S 
harmonic and decomposed into its angular harmonics. The steady flow is computed on 
acoustic mesh and provides data for the acoustic calculations. The J et > s in = lude ^^ 
steady flow potential flow model by separating the interior and extenor flow outs 'd e the 
aft fan duct with a thin barrier created by disconnecting the computationa om i • 
jet and exterior flow are allowed to merge at a defined distance downstream. In the acoustic 
Son model continuity of acoustic particle velocity is implicitly satisfied aero s the 
shear laver by careful treatment of the surface integral which appears in the finite elemen 
method '(FEM) formulation. Pressure continuity is enforced by using a penalty constraint 
on the shear layer. A model for locally reacting acoustic treatment provides a boundary 
condition on the duct walls. An attempt has been made to limit reflections on the artificial 
baffle introduced to limit the computational domain, but this is only moderately successful . 
An old but reliable frontal solution routine has been updated with considerable impact 
on computational time. Example calculations are given which show the success achieved 
in satisfying the complicated interface conditions on the shear layer and the characteristics 
of the solutions at relatively high frequencies where the refinement of the mesh becomes 
a limiting consideration for practical computations. Academic PreS s Limited 


1. INTRODUCTION 

In approach and cutback conditions the acoustic field from high by-pass ratio turbofan 
ezines is dominated by tonal noise generated by blade/vane interactions and radiated 
forward from the nacelle inlet and to the rear from the aft fan duct. In order to meet noise 
control goals active and passive techniques can be employed to control the source 
mechanisms and to attenuate acoustic propagation in the inlet and fan exhaust ducts^ 
Methods for the prediction of the effects of various noise control measures on far held 
acoustic radiation are required in the design process. The investigation reported here is 
directed toward the development of a robust computational scheme tor the prediction of 
the acoustic field attributed to tonal sources typical of blade.vane interaction in 
fan duct It is intended to be coupled to a suitable model of the source mechanism. 

The model developed is an extension of computational methods which were developed 
for inlet radiation [1-4], The inlet radiation model was based on the assumption o 


0022 ^60 X ,98/ 2202 3 5 + 23 $25.00/0/sv971480 


Q, 1998 Academic Press Limited 
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irrotational acoustic perturbations on an irrotational steady flow. A finite element code 
was developed which could accurately model the geometric details of an axisymmetric inlet 
as well as the steady flow field in and around the inlet, including the effect of forward flight. 
Rapid advances in the capabilities of work stations has made it a realistic goal to accurately 
predict the acoustic field around realistic geometries at realistic frequencies. Reported here 
is the development of a similar model for aft radiated noise. The most significant extension 
is the representation of the important effects of the fan duct jet imbedded in the 
surrounding flow which includes forward flight effects. The presence of the jet introduces 
interesting conditions which must be imposed on acoustic propagation across the shear 
layer which confines it. The methods for achieving these conditions in the context of the 
finite element method (FEM) are discussed in detail here. 

2. FORMULATION OF THE ACOUSTIC RADIATION PROBLEM 

The aft radiated acoustic field from a turbofan nacelle is described by a potential 
formulation as previously introduced for inlet acoustic radiation [1-4]. Figure 1 is a sketch 
of the important geometrical features of the aft fan duct and centre body. The nacelle has 
a forward flight Mach number M 0 . which at large distances is equivalent to a uniform flow 
directed away from the fan exhaust duct exit plane. Near the nacelle this velocity field is 
non-uniform. The exhaust flow, defined at the source plane by Mach number M,, emerges 
as a potential flow jet and extends down stream confined by a shear layer separating it 
from the exterior flow. The shear layer is terminated at a defined length at which point 
the jet and external flow merge as potential flows. The potential flow merging of the jet 
and exterior flow at the end of the shear layer produces a localized steady flow anomaly 
which has not been observed to substantially influence the acoustic radiation. 
Computations are to be carried out using the FEM in a domain including the interior of 
the aft fan exhaust duct and an exterior region made finite by invoking a radiation 
condition at an outer computational boundary and by introducing an artificial baffle 
oriented to produce a minimal effect on the acoustic radiation field. 

The nacelle geometry and the steady flow field are assumed to be axially symmetric. The 
noise source is assumed to be harmonic in time and is decomposed into its angular modal 
content, allowing a two-dimensional representation of the acoustic field in an (.v, r) plane 
through the nacelle axis of symmetry. The solution domain is shown in Figure 2. It is the 
x, r plane in cylindrical co-ordinates. The source plane is designated by C, . The fan or exit 
guide vane source is input on this plane by specifying complex amplitudes of incident duct 
modes (see references [1-4] for details of the implementation of the source boundary 




Figure 1. Sketch of the geometry of the aft fan duct, emphasizing the exhaust flow and shear layer. 
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X 


Figure 2. Computational domain showing the boundaries and regions. 


condition.). The nacelle outer surface is C„. The outer boundary of the solution domain 
Co is a circle which is a constant phase surface for an acoustic source located at the origin. 
On this boundary a radiation condition is specified. Wave envelope elements [1—4] are used 
in the far field to reach the outer boundary with minimal cost in mesh refinement. An 
artificial baffle G limits the solution domain well upstream of the fan exit plane and is 
chosen to be swept in such a way that a minimal effect on the acoustic field is created. 
This baffle is a ray from the origin and in principle at large distances from the source it 
should be non-reflecting, although near field effects do lead to reflection. The placement 
of the baffle must be considered in terms of the likely orientation of the radiated field. The 
baffle can be eliminated if computational efficiency is not a limiting factor. The shear layer 
C which separates the potential flow jet from the potential exterior flow is a rigid boundary 
for the calculation of the steady flow field and is a surface across w'hich appropriate 
continuity conditions must be satisfied in the acoustic calculations. 

The starting point for the formulation of both the steady mean flow and the acoustic 
perturbation consists of the mass and momentum equations and the energy equation in 
the form of the isentropic equation of state: 

v ■ (p\) = o, (l) 

£X + (V- V)V = -\vp, £ = (C, (2,3) 

Ot P Pn \pnj 

where p, p and V are fluid properties pressure, density and velocity, at this point in 
dimensional form, and p„ and p n are reference values of pressure and density. 
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A weak form of the field equations begins with equation (1) in which solutions for p 
and V are sought in the class of continuous functions which satisfy the weighted residual 
relation 



V W ■ (p V) - w % 
at 


dV = 


Wp \ • n dS 


(4) 


for every function W(x) in the class of continuous functions. The surface integral is over 
the boundaries of the domain of solution and n is the unit normal out of the domain. In 
the finite element discretization process which follows, the surface integral must also be 
interpreted at each subdomain boundary, namely the boundaries of the individual 
elements. The physical boundaries of the solution domain include the boundaries of the 
nacelle, including the source plane, the rigid structural boundaries, and absorbing 
boundaries such as acoustic treatment. Other boundaries are the artificial baffle introduced 
to limit the solution domain and the outer boundary of the solution domain at which a 
radiation condition is applied. All of these boundary conditions are introduced through 
the surface integral. The boundary integral is- observed to involve the mass flux normal 
to the boundary. The integral is therefore in terms of an essential conservation quantity, 
and this is typical of weak formulations in the framework of the FEM. For boundaries 
between subdomains (elements) in the FEM discretization at which there is no surface of 
discontinuity the integral produces no net contribution. This follows because on such 
boundaries the integrals on either side of the boundary produce contributions equal in 
magnitude and opposite in sign. In the present problem this applies to all boundaries 
between elements, although, as will be shown, a careful interpretation of the surface 
integral must be earned out to establish that it vanishes across the shear layer separating 
the outer flow field from the jet with a discontinuity in tangential steady flow. In particular, 
it seems to be necessary to start from the yet to be linearized form of the weighted 
continuity equation (1), and to carefully linearize it to account for the fact that on the shear 
layer, which is displaced due to acoustic perturbations, the integral is interpreted to be 
evaluated on the surface of discontinuity in tangential steady flow velocity with the unit 
normal defined to reflect this. If the field equation is linearized before the weak formulation 
is established, an essential contribution to the boundary integral is lost. 


3. BOUNDARY CONDITIONS ON THE SHEAR LAYER 


Figure 3 shows the idealized interface between the exterior flow and the jet. The surface 
of discontinuity in tangential velocity is assumed to be displaced from the mean position 
by 


0, t) = £(*, 0, t), (5) 

where £(*, 0, t) is an acoustic perturbation. The normal to the interface is tilted due to 
the slope of the shear layer approximately by 

Axih = ~j* n ' ! (6) 

where n* can be written in terms of the unit vectors n, n,, which are normal and tangent 
to the undisplaced shear layer, in the form 


n* = n 



(7) 
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Exterior tlow 


Aft fan duct 


An b 



Figure 3. Geometry of the shear layer interface, showing the acoustically displaced boundary between the jet 
and exterior flow. 


The orientation of the unit normals here are consistent with the surface below the shear 
layer for which the normal out of the fluid is in the direction of the positive normal fluid 
particle velocity, but a similar argument applies if the surface above the shear layer is 
considered. The tangential component of the normal to the shear layer is an acoustic 
perturbation quantity. 

The interface conditions across the shear layer can be determined by examining the mass 
flux and momentum flux at a moving surface of discontinuity. A particularly good 
explanation is given by Karamcheti [5]. The essential results in the case of a discontinuity 
in the velocity tangential to the surface of discontinuity are: 


(V* - V T ) • n A = (V/ - V,) ■ n* - 0, - V,) ■ n* - p,(V, - V,) • n* = 0, (8, 9) 

p u V u * n*(V. ~ V,) ■ n„ - p,V, • n*(V, - V v ) • n* = p,-p u . (10) 

Here V u and V/ are the fluid velocity above and below the discontinuity and V, is the 
velocity of an element on the surface of discontinuity. p u and p, are the fluid densities above 
and below the discontinuity and p„ and p, are the corresponding pressures. Equation (8) 
follows from the tangential component of the momentum equation and equation (9) from 
the mass continuity equation. Equation (9) is satisfied automatically due to equation (8). 
Equation (10) is the component of the momentum equation normal to the discontinuity. 
With equations (8) and (10) it is determined that 


Pu=pi , ( 11 ) 

which is the condition that pressure be continuous across the shear layer. The linearized 
version of this would require the acoustic perturbation in pressure to be continuous as well 
as the pressure of the steady flow. 

A linearized version of the surface integral of equation (4) is required for the acoustic 
analysis which follows. The fluid velocities are replaced by their perturbation forms 
V w = V r) n, 4- \ u and V/ = V ri n, 4- v, 5 where K, u and V r , are the mean flow tangential velocities 
above and below the discontinuity. The densities are replaced by p u = p ru 4- p u and 
p, = p n 4* p f with the possibly different mean densities given by p, u and p rr The acoustic 
quantities are now p u and p } . It is also important to note that the velocity of an element 
of the surface of discontinuity is an acoustic quantity and is therefore denoted by V, = v,. 
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Equation (7). and equation (8) in linearized form, are equivalent to the tamiliar conditions 
of continuity of particle displacement. 


ac , v cc 

n = t- + K„ — , 

0! CX 


, v <?c 

v, • n = — + V ,, — ■ 

Ct CX 


(12, 13) 


The linearized form of the surface integral of equation (4) on the upper and lower surfaces 
of the shear layer can then be obtained by using equation (7) and equations (12) and (la) 
(and accounting properly for the evaluation of the integral on the surtace above the shear 
layer): 



.S\, 


W(pV), ■ ru, dS = 



A\ f 


(14) 


W(p\) u * n hi dS = 


Wp r ^ dS. 
K - ct 


(15) 


It is apparent from equations (14) and (15) that along the shear layer the net contribution 
of the surface integrals will vanish if the steady flow densities above and below the shear 
layer are the same. If they are different, as in the case of a hot jet in a cold outer medium, 
there will be a net contribution which is effectively a distributed source on the shear layer 
with a strength proportional to the difference in the steady flow densities. This is 
completely consistent with equation (9). It is also consistent with the rigorous analysis 


given by Myers [6]. 




4. LINEARIZED WEAK FORMULATION 

A linearized weak formulation is obtained by continuing with equation (4) for which 
the linearization of the boundary integral has been examined. Acoustic propagation and 
radiation is modelled based on the assumption that the mean flow in and around the 
nacelle is irrotational and that the acoustic perturbation is also irrotational. The potential 
formulation makes it possible to introduce mean flow and acoustic perturbation velocity 
potentials. Acoustic perturbations are assumed on the steady mean flow such that 
$ = 4>, + <t>, p = p, + p and p = pr + p. The acoustic perturbations are assumed to be 
harmonic in time and in the angular co-ordinate such that p(x, r, d. t) = p(x, r) e""-' "" , 
p(^. r, 8 , l) * ?(*, 0 r, 8 , 0 = <p(x, r) e 1 "'' The acoustic perturbation 

in the shear layer position is also assumed to be harmonic in time and the angular 
co-ordinate yielding ftx. 9, t ) = C( v) e**' ” The steady flow density and velocity are p„ 
V<j>,. In linearized form, the weak formulation of equation (4) becomes [4] 



jjj 


\VW ■ (p r V<S) + pp>,) - i r],Wp} dV = i^ 


rr 


J J 

s, 


W( Pr , - d S 


+ 


W(p r V<p + pV(j)r) ■ n dS. (16) 


The weighting functions are taken as W(x, r , 0) = W{x, r) e' m ". Perturbations are in the 
form ofangular harmonics proportional to e"”" representing the decomposition ot the 
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solution periodic in 9 in a Fourier Series. The angular mode number is a parameter of 
the solution. The first surface integral on the right hand side is on the shear layer S K and 
the second surface integral is over all remaining surfaces bounding the domain. Notice that 
the unit normal for the second integral is the normal out of the domain at the surface in 
question. The weak formulation continues with the linearized momentum equation [4] 

p = — p (i t] r (j) + V(p r ■ V(p), (17) 

which is used to replace p in equation (16), the linearized equation of state, 

p = c;p, ( 18 ) 

which is used to produce an alternative form of the momentum equation in terms of 
acoustic pressure, 

p = — p r {\r] r (j) + V<f) r ■ Vcp). (19) 

Equation (19) is used to define acoustic pressure difference on the shear layer and to 
post-process the field solution for (p to obtain the acoustic pressure field. The acoustic 
particle velocity and acoustic velocity potential are related according to 

v = Vcp. (20) 

The linearization process also produces the weighted residual formulation for the steady 
flow, 


VW • (p r V<j)r) dV = 


W(p r V(p r ) * ndS, 


and the steady flow momentum equation in terms of the speed of sound, 

c; = 1 - [P0, • V<j) r - M{], 

and in terms of the steady flow density, 


Pr = 



(7-1) 


“ 1 ly - I ) 

(V(pr * V<P r ~ 


( 21 ) 


( 22 ) 


(23) 


Equations (16) through (23) are in non-dimensional form where (p is the acoustic potential, 
(p r is the local mean flow (reference) potential, p is the acoustic density, p r is the local mean 
flow density, and c r is the local speed of sound in the mean flow. All quantities are made 
non-dimensional by using the density in the far field, the speed of sound in the far 
field, c % , and a reference length which is defined as the duct radius at the source plane. 
/?, where acoustic modal amplitudes are defined. This plane could be the fan plane or the 
exit guide vane plane, but it is not restricted to these locations. The acoustic potential is 
non-dimensional with respect to c*R, and the acoustic pressure with respect to p,.c\. 
Lengths are made non-dimensional with respect to R. Time is scaled with Rjc^, leading 
to the definition of non-dimensional frequency rj r = (vRlc\; oj is the dimensional source 
frequency and ;V/ X = M„ is the Mach number in the far field representing the forward flight 
effect. 
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Equation (21) is the weighted residual formulation for the calculation of the 
compressible potential flow within and around the nacelle. Equations (22) and (23) are 
subsidiary relations that can be used in an iterative solution which at each stage uses a 
density field derived from the previous iteration step. F0 f , c r , p r are required data for the 
weighted residual formulation of the acoustic problem. In the results reported here only 
the first iteration of this process is used to define the potential flow field. This is 
accomplished by solving the incompressible problem and then computing a variation in 
steady flow density and speed of sound. 

The second surface integral in equation (16) provides the boundary conditions on the 
duct walls and on the source plane. The modelling of duct acoustic treatment in the context 
of this integral is discussed in a later section. The acoustic source is specified by the complex 
amplitudes of acoustic duct modes at the source plane. On this plane the FEM modal 
values of acoustic potential are replaced by the complex amplitudes of the acoustic 
potential modes by an eigenfunction expansion. The incident acoustic modal amplitudes 
are input and the reflected modal amplitudes are computed as part of the solution. Details 
of this procedure are available in references [1-4]. 

The same surface integral provides the mechanism for introducing the boundary 
conditions on the artificial baffle and the non-reflecting boundary condition at large 
distances on the outer boundary of the computational domain. These details are also 
explained extensively in references [1-4]. 

Acoustic pressure and particle displacement are continuous across the shear layer. The 
continuity of particle displacement is implicit in the handling of the surface integral on 
the shear layer. Continuity of pressure must be explicitly enforced. The implementation 
of this condition will be discussed presently. 


5. COMPUTATIONAL MESH 

A particularly sensitive issue which must be resolved is the construction of a mesh which 
is consistent w’ith the geometry requirements and which can be generated automatically 
from data describing the nacelle and centre body. It is essential that the mesh be structured 
to minimize the bandwidth for the linear equation solver. The major constraining feature 
is that the trailing edge of the fan duct is thin or cusped. In addition, the near field mesh 
must evolve into a smooth transition to the far field wave envelope mesh. 

In order to meet all of these requirements, a mesh which combines quadrilateral and 
triangular elements has been used. Figure 4 shows the details of the near field mesh. The 
interior of the duct and the extended jet uses conventional eight-node quadrilateral 
elements. Most of the exterior region is also composed of quadrilateral elements. However, 
a fan shaped region of six-node triangular elements is used to allow a transition around 
the sharp trialing edge. Primarily due to the constraint on the mesh structure for 
minimizing bandwidth, this transition would not be possible with rectangular elements 
without introducing severe distortion in the neighbourhood of the trailing edge. The far 
field mesh utilizes the wave envelope element concept [1-4], and presents no problems. A 
relatively coarse near and far field mesh is shown in Figure 5 and the wave envelope 
element region can be seen. Note in Figure 4 that the exhaust duct trailing edge is reflexed, 
representing the most severe case. 

Mesh generation produces two mesh connectivities. For the potential flow code, velocity 
potential is discontinuous across the shear layer dividing the extended jet from the external 
flow. Elements are therefore disconnected across the shear layer. For the radiation code, 
acoustic velocity potential is also discontinuous across the shear layer. Elements above and 
below the shear layer have additional degrees of freedom on the shear layer boundaries 
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representing acoustic particle displacement, which is continuous across the shear layer. The 
mesh for the radiation code therefore introduces two degrees of freedom at the nodes on 
the shear layer. There are 1 1 degrees of freedom for rectangular elements and nine degrees 
of freedom for triangular elements on the shear layer. Pressure continuity is enforced by 



Figure 5. Aft fan duct far field mesh with wave envelope region. 
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Figure 6. Bridging elements on the shear layer. 


using a penalty constraint on the shear layer, and it is not necessary to introduce pressure 
as an additional variable on the shear layer. However, it has been found convenient to 
introduce six-node “rectangular” transition elements ot zero thickness between the 
standard elements above and below the shear layer for generating the “penalty element 
stiffness matrices”. In order to maintain consistency in the meshes for the potential flow 
calculations and radiation calculations, these elements are accounted for in mesh 
generation for both codes. Details of the elements on the shear layer are shown in Figure 6 
where for an example a bridging element is inserted between a triangular element above 
the shear layer and a rectangular element below the shear layer. 


6. STEADY FLOW CALCULATIONS 

A potential flow code generates the steady flow field in and around the aft fan duct. 
Incompressible potential flow has been assumed as a first approximation. Variations in 
density and speed of sound are based on the isentropic equation of state and 
incompressible velocity field with specified conditions on Mach number, density, and speed 
of sound in the far field. It is within the framework of the general formulation to treat 
the potential mean flow as compressible, and the present approach can be viewed as the 
“zeroth” iteration of the compressible isentropic case. Fully compressible isentropic mean 
flow has been used by the first author in acoustic propagation and scattering calculations 
in pipes in which no far field radiation is modelled. In the type of problem considered here 
the computational overhead required for the several iterations necessary to produce a 
compressible mean flow has not been considered justifiable at the present stage of 
development. 
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Figure 7. Potential field for the steady flow from the aft fan duct and in the surrounding flow field. This case 
corresponds to M, = 0-5 and forward flight velocity \f„ = 0-2. 


The potential flow field has been structured to include flow in a jet region downstream 
of the fan duct exit plane. This has been done with the introduction of a '"rigid” duct 
boundary representing the fan exhaust shear layer which extends the prescribed length of 
the jet. The ngid boundary is introduced by permitting the velocity potential to be 
discontinuous across the shear layer. At the downstream end of the shear layer the 
discontinuity in velocity potential is terminated and merging of the interior and external 
flows is permitted. The merging can produce very high velocities and reverse flow near the 
termination of the shear layer. This is smoothed out by restricting the velocity near the 
end of the shear layer to neither go above a reference velocity which is determined midway 
along the underside of the shear layer (in the jet) nor to go below a similarly determined 
velocity on the upper side of the shear layer (in the outer flow). In the near field mesh of 
Figure 4 the shear layer boundary can be seen to extend downstream about two duct radii. 
The merging distance is adjustable, and is chosen to provide sufficient distance for full 
effect on the acoustic radiation, and to move the perhaps unrealistic merging region away 
from the important part of the acoustic field. Computations for the steady flow field are 
carried out on the same mesh used in the acoustic case. This is done so that the steady 
flow data is produced in a form compatible with the acoustic mesh. The mesh is invariably 
much more dense than would be required for the steady flow calculations, however, the 
problem is symmetric, and the solution routine is considerably faster than for the 
comparable acoustic problem (about three or four times faster for large meshes). 

A typical potential flow field is shown in Figure 7, where the jet and surrounding forward 
flight effect contours of velocity potential are clearly shown. 


7. ACOUSTIC PRESSURE CONTINUITY ON THE SHEAR LAYER 

The linearized weak formulation of equation (16) has the subsidiary condition of 
continuity of acoustic pressure across the shear layer. This condition is not easily satisfied 
because the formulation is in terms of acoustic velocity potential. However, equation (19) 
provides a connection between the acoustic velocity potential above and below the shear 
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layer which can be exploited to implement the continuity condition by using a penalty 
method [7]. 

The important features of the penalty method can be described relatively easily. 
Equation (19) and the condition pressure continuity on the shear layer are used to obtain 


Ap = p„ -p, = pr\ i n,(p, + Mi 


— pA \r\ r (p„ + M „ 


caK 

5x 


0 . 


(24) 


Equation (24) applies computationally on the shear layer, r = r,, x { ^ x ^ -\r : , where r, is 
the radius of the axially symmetric shear layer and X\, xi are the axial co-ordinates ol the 
left and right ends of the shear layer; jti coincides with the trailing edge of the nacelle at 
the fan exit plane. The subscripts / and u denote values ot the steady state and acoustic 
quantities below and above the shear layer. In the finite element context, the acoustic 
potentials and (j> u can be written notionally in terms of a global interpolation matrix 
(yV(.x)]. For example, 


4>'(x) = ( 25 ) 

where (j) f is the vector of nodal values of (pi(x) below the shear layer. The interpolation 
matrix [iV(,x)] is a row matrix with elements N,{x), / = 1 , AW, where AW is the number 
of finite element nodes and N,(x ; ) = 1, / =y, Ni(Xj) = 0, / ^ j, j — 1, AW. The substantial 
derivative operators in equation (24) are defined such that, for example. 


) = pr,[ i tfr -h Ml 


cx 


4 »- 


(26) 


In vector-matrix format. 


D,(<l> t ) = [N{x)][D t ]<fi t , (27) 

where [D f ] is a diagonal matrix of operators p r/ (nF + Mi dldx). Equation (24) can be written 
as 


M(*)][0f. M T = mx)][Di] - [N(x)][D u ]M>, 4>»Y = 0. (28) 

The modified interpolation matrix [W^x)] has elements which can be viewed as 
interpolating Ap(x) from nodal values of the acoustic potential on either side of the shear 
layer. The weighted residual of equation (28) is formed on the shear layer using as 
weighting functions elements of [N* (,x)], which are the complex conjugates ol the 
interpolation functions. This yields 


f [N?(x)mw] d S,[<f>i, <t>A T = 0. (29) 

J J 

s, 

This is a weighted residual (or variational) statement that Ap vanishes on the shear layer. 
It produces a “stiffness matrix” which is consistent with this statement. If this is appended 
to the weighted residual formulation of equation (16) with a large multiplier a, it forces 
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the weighted residual formulation to have a solution constrained by equation (24). The 
modified weighted residual statement is written as 


[VW‘(p r V<p + pV(j> r ) - i rj,Wp} dV = U]r 


W{ Pn - pJC d5 


+ 


W(p r V (\> + pV(j)r) • n dS 


Wdp dS. 


(30) 


The weighting functions W{x) introduced in equation (30) are the pressure difference 
interpolation functions identified in equation (28). The penalty integral, equation (29), is 
introduced along the shear layer and produces penalty stiffness matrices which bridge the 
shear layer and include nodal values of acoustic potential on both sides ot the shear layer. 
This is most easily implemented in the finite element context by introducing transition or 
bridging elements on the shear layer as shown in Figure 6. These elements are ot zero 
thickness with three nodes on the top and three on the bottom to connect to the three nodes 
on the conventional triangular or quadrilateral elements above and below the shear layer. w 
Finite element assembly proceeds as with other elements in the mesh. No new global nodes 
are introduced and there is no change in the bandwidth of the formulation nor to the 
sequence of operations in the equation solving procedure. 

The boundary integral on 5 represents natural boundary conditions which must be 
imposed on the other boundaries of the domain. The far field boundary C x is at a large 
distance from the nacelle and is a non-reflecting surface on which a radiation condition 
is applied via the boundary integral. This surface is the outer boundary of wave envelope 
elements which allow a transition from a fine mesh near the nacelle to a very coarse mesh 
in the far field. Most of the nacelle and centre body surfaces are rigid, where the normal 
component of acoustic particle velocity vanishes. In addition, the normal component of 
the mean flow velocity also vanishes and the flow tangency condition requires that 
Vcfrr - n = 0, eliminating the boundary integral. A portion of the fan duct and centre body 
is acoustically treated. On these surfaces an impedance relation is specified, and this can 
be introduced through the boundary integral. The acoustic source is also introduced using 
the boundary integral. Details of the imposition of natural boundary conditions can be 
found in references [1-4]. 

In the results presented here there is no difference in steady flow density across the shear 
layer. The boundary integral on S* which arises from considerations of conservation of 
acoustic particle displacement across the shear layer therefore vanishes. 


8. ACOUSTIC TREATMENT ON DUCT WALLS 

In the FEM formulation described here provision has been made for acoustic treatment 
on the duct wall and on the centre body. In the present context the acoustic field is 
described in terms of an acoustic potential formulation, while the boundary condition 
relates pressure and particle velocity. The implementation is described in this section. 

A locally reacting acoustic lining material specified by its frequency dependent 
impedance or admittance is placed on an interior surface of the duct. The boundary 
integral of equation (30) is the mechanism by which the boundary condition imposed by 
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this acoustic treatment is introduced. On surfaces where acoustic treatment is present the 
normal component of mean flow velocity vanishes and the lining boundary integral 
simplifies to 


Il = 


Wp,V(j> • n d S, 


N. 


(31) 


where v • n = V<\> • n is the normal component of acoustic particle velocity, v. The unit 
normal n is directed out of the computational domain and therefore into the acoustic 
treatment. The acoustic treatment is described by a local impedance relationship which 
connects acoustic pressure to a conceptual wall displacement velocity. In general, the types 
of acoustic treatment of interest are porous and the wall intself docs not displace but the 
fluid in the pores does. It is the fluid velocity in the porous wall, directed normal to the 
wall, which is referred to as wall displacement velocity. The impedance relationship is of 
the form 

£ = Z = T (32) 

V n A 

where p is the non-dimensional acoustic pressure and v„ is the non-dimensional wall 
displacement velocity, directed into the wall. The impedance Z is a prescribed function of 
frequency and is non-dimensional with respect to p x t*x, that is, the dimensional impedance 
would be p x c x Z . A is defined as the non-dimensional acoustic admittance. The relation 
between the fluid particle velocity at the wall and the wall velocity is one of continuity 
of particle displacement. This yields 


v 


n = 4- 


at ox 


(33) 


where £(x, 0, /) = £(.v) e ' 1 "' l ~ m0) is the wail displacement normal to the wall mean position, 
positive directed into the wall, and related to v„ by v n = A/„ is the steady flow velocity 

at the wall, non-dimensional with respect to . The choice of positive or negative sign 
depends on whether the acoustic treatment is on the outer or inner wall of an annular duct. 
It is assumed that all lined surfaces have negligible curvature in the direction of the duct 
axis so that the rigorous description of the flow/surface kinematics [6] is simplified. With 
harmonic time dependence. 


v • n = -r- \ \rj r 4- M w ]i\ t . 
I YJ r \ ox f 


(34) 


The relation between acoustic particle velocity and acoustic pressure is 


v n = ±( \T} r + ^ )Ap. 


(35) 


The relation between acoustic pressure and acoustic velocity potential is provided by the 
acoustic Bernoulli equation of equation (19). Equation (35) can be rewritten 


v n = + — [ irjr -b M w 7 - 
lrjr \ OX 


p r A[ ir! r + 


(36) 
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The boundary integral becomes 


4- prWv * n dS = i rj r 


WAp;<p dS + 


WAp;M„ ^ dS 

OX 


+ 


C.v rjr 


Wp r Mri~ 

cx 


Ap,M,M 

ox 


dS. 


(37) 


The upper and lower sign choice depends on whether the outer or inner wall integral is 
considered. The first two integrals on a boundary where acoustic treatment is present are 
easy to implement in the finite element formulation because only continuity of acoustic 
potential is required. The admittance, A , is assumed piecewise continuous and non-zero 
on a portion of the interior surface of the duct. An integration by parts, which is essentially 
an application of Stokes’ Theorem on the interior surface, is performed to make the last 
two integrals compatible with the weak formulation. The integral representing the 
boundary condition on interior surfaces can now be written 


+ 


p r W's • n dS = i Y\ f 


WAp;<f) dS -h 


WAp;M v dS 
ox 


Sl 


(Ap4)-^{Wp,M„)dS + - 


A p r M w 


d<fi 


ox 


( Wp r M w ) dS. (38) 




Equation (38) is in a form which is appropriate for application of standard finite element 
techniques to generate “boundary matrices” which are appended to the element stiffness 
matrices of elements whose outer boundaries represent acoustically treated surfaces. 


9. AN ABSORBING BAFFLE 

A restriction of the present FEM mesh is the presence of the baffle which is used to limit 
the computational domain, presumably with little reduction in the quality of the solution. 
It is assumed that the baffle is swept back at least 90 c from the angle of peak radiation, 
however, this condition is often violated because it requires a mesh generation change to 
accomplish it. In theory the baffle is non-reflecting at large distances from the nacelle since 
it is a ray extending from the origin [2]. Near the intersection of the baffle and the nacelle 
the baffle is a reflecting surface and its presence has the possibility of contaminating the 
solution with spurious reflections. Experience has shown that the baffle has little effect on 
the peak lobe in the radiation pattern if the 90° rule is adhered to. Flowever, there has 
been interest in using the inlet code and aft radiation code to generate the SPL directivity 
on the full 180° arc around the engine. This would be accomplished by separately obtaining 
the inlet and aft radiation results and then superposing them. Presumably, the peak lobes 
fore and aft would be little affected but the region at 90° to the engine axis would be 
critically dependent on a legitimate superposition of the inlet and aft results. This is not 
possible to achieve because of the baffle, unless it is completely non-reflecting. 

An investigation has been made of the possibility of making the baffle at least partially 
non-reflecting. This has been done by introducing absorption on the baffle. As in the case 
of the nacelle acoustic treatment, this is done through the surface integral on S in 
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equation (30). On the baffle it is assumed that the flow is adequately represented by the 
uniform Mach number = M*i. This is true far from the nacelle, and is approximately 
true near the nacelle. The acoustic density perturbation is given by equation (17) evaluated 
with p r = 1 and c, = 1, assuming that far field steady flow conditions apply on the baffle. 
The surface integral on the baffle can then be written as 


W(p r V(f) + pP(j>r) ' n dS = 


w 


V<p - n — (M x * n) 


d 4 + M, 

ct 


V(p 


d5, (39) 


where n is the unit normal of the computational domain. The impedance condition on the 
baffle surface is defined simply as 


p = — v • n. 

Z*/p,c« is the non-dimensional impedance and Z h jp*c^ = 1 iA, where A is the 
non-dimensional admittance. This impedance condition corresponds to no real physical 
situation but rather is introduced to provide absorption on a notional boundary through 
which there is a steady mean flow. The acoustic Bernoulli equation (19) and the definition 
of the acoustic velocity potential 

v = V(p ( 41 ) 

leads to the conclusion that on the baffle, 


V(j> • n = 


gx£x 

Z* 


(f +M -- p 4 


(42) 


The boundary integral can therefore be written as 


W(p r V(j> + pP<M • n d5 = — 


+ n )( ^ + M x • V<p 


dS. 


(43) 


The boundary integral of equation (43) is applied only in the near field portion ot the baffle. 
In the wave envelope region the theoretically non-reflecting character of the far field baffle 
is left unchanged. The introduction of a locally reacting impedance boundary on the baffle 
cannot be expected to produce complete absorption any more than on the wall ot a duct. 
As will be shown, only a modest absorption can be achieved. 


10. POSTPROCESSING 

Postprocessing of the acoustic velocity potential solution using the acoustic Bernoulli 
equation (19) to obtain acoustic pressure can be carried out in two ways. The approach 
which is most efficient computes acoustic pressure at the element nodes using the element 
shape functions. The nodal values are then averaged, to account for the fact that 
derivatives of potential are not continuous across element boundaries in the FEM 
formulation. For sufficiently fine meshes this produces acceptable results. Results presented 
in this paper are obtained by this method. 

A second method available for postprocessing acoustic velocity potential to obtain 
acoustic pressure carries out the calculations at Gauss points in each element. The Gauss 
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points are known to be points at which optimal accuracy is achieved in the calculation 
of derivatives and therefore in calculation of acoustic pressure. The number of Gauss 
points is generally less than the number used in the Gauss integration in the formulation 
of the element stiffness matrices. The array of solution points on the grid constructed in 
this way can then be plotted as contours of equal acoustic pressure or equal sound pressure 
level using one of several available commercial plotting packages. 


11. SOLUTION TECHNIQUES 

The principal advantage of the FEM formulation described here is that it is 
computationally relatively efficient and therefore provides a useful tool for design 
calculations. This efficiency decreases as the non-dimensional frequency of the acoustic 
source increases, requiring a proportional increase in the mesh density and a 
disproportionate increase in computation time (by approximately the square of the ratio 
in mesh density). For this reason it is appropriate to give some observations on the linear 
equation solving routine which accounts for almost the entire computational time. 

Previous fan noise radiation codes [1-4] used a frontal solver due to Irons [8]. This was 
extremely efficient in the use of active memory, partly because of considerable data transfer 
using direct access I/O in storing and retrieving element stiffness matrices in the 
assembly/solution procedure and in retrieving mesh and steady flow data. The resulting 
direct access files were subsequently read many times in the various FEM operations and 
in postprocessing. This efficiency in storage was offset by a significant cost in execution 
time. Nacelle design and source modelling have become the primary uses of the codes and 
execution time is a primary issue in a work station environment in which storage has 
become a much less limiting factor. Direct access operations are efficient from a 
programming standpoint, but inefficient in I/O time. In the version of the radiation code 
reported here all direct access I/O has been eliminated in favour of active storage or 
sequential I/O. This has resulted in as much as 50% reduction in computation time, 
dependent mainly on available fast memory. 

Experiments with several popular iterative solution routines show that for the 
two-dimensional structure of problems considered here the direct solvers are always faster. 
This is consistent with the experience of other investigators [9], There are indications that 
iterative solvers can be faster for similar problems in a three-dimensional geometry. The 
choice has been made to retain the modified Irons frontal solver. 


12. EXAMPLE CALCULATIONS 

The example calculations shown here are obtained on a mesh with about 27 000 degrees 
of freedom which is shown in Figures 4 and 5. This mesh becomes inadequate for 
non-dimensional frequencies much in excess of rj r — 25, and with the element distribution 
shown does better for acoustic radiation toward the sideline (high angular mode number 
or high radial mode number). Angular mode order corresponds to the value of m in the 
angular Fourier component e - ""* 1 . Radial modes for a specified angular mode are 
enumerated beginning with n = 1. The geometry of the aft fan duct is generic, including 
an extended centre body and thin fan duct lip, in this case reflexed. The exterior Mach 
number is M„ = 0-2 and the jet Mach number at the source plane is M, = 0-5. 

The first result shows the success of the penalty method in implementing the condition 
of continuity of acoustic pressure across the shear layer. This is most effectively shown at 
low frequency because the acoustic field is relatively simple and the discontinuity in 
acoustic potential and continuity in acoustic pressure is easy to see. The frequency chosen 
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Figure 8. Contours of equal acoustic potential with \f, = 0*5, Mo = 0-2. Reduced frequency t] r — 5*0, input 
angular mode m - 2, radial mode ^l,no acoustic treatment. Acoustic potential is discontinuous across the 
shear layer. 


is rj, = 5 with an angular mode m = 2 and radial mode n = 1 input with unit pressure 
amplitude. The mesh is quite adequate for this low frequency. Figure 8 shows the radiated 
field in terms of contours of constant acoustic potential magnitude superimposed on the 
computational domain. In this example only five contours are produced to simplify the 
plot. The contours range from 15 dB above the maximum level on the outer boundary to 
15 dB below. In Figure 8 it is clearly seen that the acoustic potential is discontinuous across 
the shear layer. Figure 9 shows similar contours of acoustic pressure and these are seen 
to be continuous across the shear layer. The pressure has been obtained by post-processing 
the potential field by using equation (19) with FEM interpolation at the nodes. Pressures 
thus obtained are averaged at common nodes. It is important to note that nodes across 
the shear layer are not common and the pressure across the shear layer is not averaged. 
The effectiveness of the penalty method is demonstrated by this example, as is the quality 
of the solution at this low frequency. Figure 10 shows an additional method of presentation 
of the radiation directivity. This represents calculations of sound pressure level on a 
circular arc at a radius of 10 duct radii from the origin, normalized to 100 dB maximum. 
In this case it emphasizes how broad the principle lobe is near the peak. 

The results of Figure 10 can be used to compare the peak radiation angle in the principal 
lobe in the far field to predictions obtained using ray theory. A code has been written which 
is based on the analysis of Rice and Saule [10] for estimation of the radiation directivity 
from an exhaust duct. It is based on an extended analysis since it considers annular ducts 
while the original work of Rice and Saule considered only circular ducts. Propagation 
angles in the duct are determined from a formal eigenvalue/eigenfunction analysis and the 
convection and refraction effects are included as in reference [10]. It is predicted that the 
group velocity in the duct at the specified frequency and in the specified mode, rj, = 5, 
m = 2, n = E is at 34-7°, and the phase velocity is at 51.2°. The peak propagation angle 
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in the far field is estimated to be 49*5\ The observed angle of peak radiation in Figures 
7 and 8 is around 50°, but the peak is so broad that it is difficult to pick the angle precisely. 
The correlation is excellent, although it must be pointed out that the mesh origin (0 5 duct 
radii back from the duct exit plane) is used in defining the directivity in this example. The 
Rice/Saule analysis would be based on an origin at the duct exit plane. Because the peak 
lobe is so broad there is little point in examining the effect of the origin shilt on the stated 
comparison. This will be done in the next example which produces a sharper peak lobe. 

A higher frequency case with a lower angle ot peak radiation is the second example. 
In this case the non-dimensional frequency is r\ f = 25 and the modal input is m — 10, n = 1 . 
This is getting close to the limit of resolution for the mesh. Figures 11 and 12 show the 
two types of presentation for acoustic pressure. Figure 1 1 showing contours ol constant 
SPL, while generally reasonably clean, emphasizes the assertion that the limit ot resolution 
is close at hand. The breakdown of the mesh adequacy always appears in the near to 
intermediate field first and is usually related to mesh density in the region between the near 
field and the wave envelope region. The number ot elements required in the generally radial 
direction is critical, and this can be minimized by bringing the wave envelope region in 
as close as possible. In the aft radiation case the jet interferes with this, and the wave 
envelope region must start far enough out that the jet is nearly entirely merged with the 
exterior flow. Figure 12 shows the polar directivity based on an origin at the exit plane 
(non-dimensional x = 0-5) and demonstrates that these far field calculations are generally 
better than the near field because of the wave envelope interpolation. This mesh has been 
pushed to Y\ r = 35 without complete failure, and has the characteristic that it produces 
better results for modes which radiate well to the sideline than for those which radiate at 
relatively small angles to the axis as in these examples. This probably results from the 
complicated interaction of transmission and reflection of modes at near grazing incidence 



Figure 1 1. Contours of equal acoustic pressure with M, — 0-5, Afo = 0-2. Reduced frequency r\ , -5 0, input 

angular mode m = 10, radial mode n = 1, no acoustic treatment. Acoustic pressure is continuous across the shear 

layer. 
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at the shear layer. The mesh shown in Figures 4 and 5 has proven to be a good generic 
structure to work with. 

Figure 12 can be used to compare the peak radiation angle in the principal lobe in the 
far field to predictions obtained using ray theory. It is predicted that the group velocity 
in the duct at the specified frequency and in the specified mode, rj r = 25, m = 10, n = 1, 
is at 27*4 C , and the phase velocity is at 39-4°. The peak propagation angle in the far field 
is estimated to be 43*9\ The observed angle of peak radiation in Figure 12 is around 42' 
and is adjusted for the origin shift to the exit plane. The correlation with the Rice Saule 
result is good, particularly when it is noted that flow conditions along the shear layer in 
the FEM calculations are not uniform, and within the jet region the Mach number is 
reduced below M, = 0-5 due to the gradual reduction in radius of the centre body. The 
effect can be observed if an average Mach number M, = 0*45 in the jet is used in the 
Rice/Saule approximation. The ray prediction would yield 4l*9~ which is about the same 
as the FEM prediction which accounts for the non-uniform Mach number in the jet. 

An example of the effect of acoustic treatment on the duct walls is also shown in Figure 
12. A locally reacting lining with normalized impedance and admittance given by 
Z = 1*221 - iO- 122, A = 0*81 1 -f i0*08 1 is assumed in the high frequency case. The 
impedance/admittance is optimum for the rj r = 25, M, = 0*5, m — 10, n = l mode. 
The outer wall of the fan exhaust duct and the centre body are lined over a length oi 
0*916 R beginning at 0*074 R forward of the assumed source plane. The attenuated 
directivity shown in Figure 12 reveals an attenuation of as much as 5 or 6 dB at polar 
angles below the shifted principal lobe which is now at about 45\ What was once a 
relatively well-defined principal lobe is now considerably broadened and beyond 45' there 
are areas in which the SPL is increased, primarily due to filling in of interference notches. 
The angle shift of the principal lobe is consistent with the fact that the effect of the acoustic 
treatment would be to increase the angle of the phase velocity and group velocity vectors 
(lower the cut-off ratio) within the duct. 

Finally, Figure 13 is used to show the effect of an attempt to reduce the effects of 
reflection from the baffle. A resistive ‘dining" with non-dimensional admittance 
A = 0*8 4- iO-O has been placed on the baffle in the region of conventional elements (the 
wave envelope elements in theory should not produce reflections). It is seen that there is 
an observable change in SPL at large polar angles (near the baffle) and in the region near 
the exhaust axis where the directly radiated field is at low SPL. Since it is not known what 
the true reflection free directivity should look like, no conclusion can be drawn other than 
the baffle does have an effect on the directivity at large angles, and that the effect is 
modestly changed when the baffle is made dissipative. Perhaps of more importance is the 
fact that virtually no effect is observed near the principal lobe, suggesting that the baffle 
has correctly been assumed to be non-intrusive in this region. 

While not entirely definitive, the results shown here suggest that the FEM model of aft 
fan radiation captures the known refraction effects of the shear layer very well. Extensive 
bench marking of the code against experiments has been carried out by other investigators 
[11, 12]. Comparisons of calculations and measurements have been very good. 


13. CONCLUSION 

A finite element model for acoustic propagation and radiation within and exterior to 
the aft fan duct of a high by-pass turbofan engine has been developed. It is based on the 
assumption of irrotational acoustic perturbations on an irrotational steady flow. The jet 
is modelled in the steady flow calculations by a potential flow constrained by a shear layer 
and allowed to merge with the surrounding flow downstream of the fan duct exit plane. 
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The formulation is restricted to axisymmetric geometries and harmonic sources described 
by their angular and radial modal content. The condition of acoustic particle displacement 
continuity across the shear layer is shown to be satisfied by proper interpretation of a 
boundary integral which occurs in the FEM formulation. Continuity ot acoustic pressure 
is implemented by introducing a penalty method based on the relationship between 
acoustic pressure and velocity potential. Example calculations have shown that the 
continuity of pressure is accurately enforced. Resolution ot accurate solutions at high 
non-dimensional frequencies is limited by mesh density which has implications on storage 
requirements and execution time. In the present study computations with over 27 000 
degrees of freedom have been shown to produce reasonable results up to the reduced 
frequency rj r = 25. Doubling the frequency would require an approximate doubling of the 
density of the mesh. 
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Variable order mapped infinite wave envelope elements are developed for 
finite-element modelling (FEM) of acoustic radiation in a uniformly moving 
medium. These elements can be used as a non-reflecting boundary condition for 
computations on an infinite domain in which a radiating body is immersed in 
a moving medium which is essentially undisturbed outside of the near field. An 
additional result of this study shows that the mapped wave envelope elements 
provide a boundary condition equivalent to stiffness, mass, and damping matrices 
appended to the inner mesh. By choosing the transition between the standard FEM 
mesh and the mapped infinite wave envelope as a surface of constant phase the 
mass matrix is caused to vanish identically. This has implications for transient 
FEM modelling of acoustic radiation. A demonstration of the characteristics of 
mapped infinite wave envelope elements is given in the context of acoustic 
radiation from a turbofan inlet for which benchmark results are known, 

C 1999 Academic Press. 


1. INTRODUCTION 


Modelling of acoustic radiation is usually complicated by the requirement that 
prediction of the acoustic field is required in some finite subdomain of an infinite 
domain. This requires that computations be limited to the subdomain with 
a non-reflecting boundary or that the infinite domain be mapped on to a finite 
computational domain. In finite element modelling this has led to the study of 
a number of forms of infinite elements [1-3], wave envelope elements [4, 5], and 
mapped infinite wave envelope elements [6-9]. The several forms of infinite 
elements in some sense map the infinite domain to a finite domain. Wave envelope 
elements restrict computations to a large but finite domain bounded by 
a Sommerfeld radiation condition. The non-reflecting boundary is reached from an 
inner standard finite-element domain via large elements in which the shape functions 
are augmented to reflect decay with distance from the source and the temporal and 
spatial character of outgoing waves. The attnbutes of infinite elements and wave 
envelope elements are combined in mapped infinite wave envelope elements. 

0022-460X/99/290665 + 23 S30.00/0 C 1999 Academic Press 
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Mapped infinite wave envelope elements have been investigated extensively for 
acoustic radiation in a stationary medium. They have certain apparent advantages 
as compared to standard wave envelope elements. In the case of harmonic 
radiation, the most significant advantage is the possibility of adjusting the order of 
the elements to fit the requirements of the problem. Formulation of the elements 
reveals the possibility of including within the element shape function an explicit 
dependence on inverse powers of the distance from the apparent acoustic source, 
consistent with theoretical results [8], This fact allows the introduction of mapped 
infinite wave envelope elements well into what would normally be considered the 
acoustic near field, reducing mesh refinement and dimensionality. The shape 
functions in mapped infinite wave envelope elements can accommodate nearfield 
effects, and this fact can be enhanced by adjusting the order of the interpolation in 
the elements to fit the problem requirements. A second advantage in the FEM 
formulation using the mapped infinite wave envelope elements is that the 
Sommerfeld radiation boundary is infinitely far away and is never explicitly 
appended as a natural boundary condition. Astley et al. [9] also demonstrated the 
applicability of mapped infinite wave envelope elements to problems in transient 
acoustic radiation, a feature which has not been exploited in standard wave 
envelope elements. With an appropriate choice of mesh geometry they show that 
mapped infinite wave envelope elements provide a boundary condition which is 
well suited for time-marching solutions. The advantages of the mapped infinite 
wave envelope elements are not without cost, and the trade-off comes in the form of 
increased band width of the discretized field equations that is introduced by high 
order mapped elements. This may offset efficiency gains achieved by reduction of 
the extent of the computational near field and therefore the standard FEM mesh if 
bandwidth-sensitive solvers are used. 

The study reported here extends the variable order mapped infinite wave 
envelope concept to uniform steady flows, principally in connection with 
aeroacoustic problems related to turbofan acoustic radiation. This is a direct 
extension of the development of Astley et al. [8, 9], They present their formulation 
in the context of problems in three dimensions in Cartesian co-ordinates. The 
application here is in a cylindrical co-ordinate system reduced to two dimensions 
by taking advantage of periodicity of the solution in the angular co-ordinate. The 
development of the mapped wave envelope elements is completely general and not 
restricted to this co-ordinate system. Harmonic radiation is considered explicitly; 
however it is shown here that as in the case of a stationary medium, with a judicious 
choice of the mesh geometry, the structure of the mapped elements becomes 
favorable for transient calculations. 


2. AN APPLICATION TO TURBOFAN INLET ACOUSTICS 

An important problem of acoustic radiation in a moving medium is available in 
the study of the acoustic field of a turbofan inlet. The noise due to turbo-machinery 
sources within the inlet is propagated in the inlet and radiated to the (infinite) far 
field. Acoustic propagation and radiation occurs in a high-speed potential flow 
which is the net effect of flow into the inlet and the forward flight of the inlet. In the 
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r 



Figure 1. Computational domain showing genetic geometry of 
regions used in the finite element, wave envelope, and mapped 
formulations. 


the nacelle and boundaries and 
infinite wave envelope element 


steady flow far field (perhaps nearer to the inlet than the acoustic far field) acoustic 
radiation occurs in a uniformly moving medium. It is required to make 
computations to predict the radiated field in a finite subdomain relatively near the 
inlet. This has been approached in the past by terminating the computational 
domain with a Sommerfeld condition on a boundary reached by the use of wave 
envelope elements [5, 10-12]. Here it is intended to investigate the application of 
mapped infinite wave envelope elements to obtain closure of the computational 

domain. 

For turbofan inlet acoustic radiation the nacelle geometry and the steady flow 
field (representing flow into the inlet and forward flight) are assumed to be axially 
symmetric. The noise source is assumed to be harmonic in time and is decomposed 
into its angular modal content, allowing a two-dimensional representation of the 
acoustic field in a plane through the nacelle axis of symmetry. The solution domain 
is shown in Figure 1. It is the x, r plane in cylindrical co-ordinates. The source 
plane is designated bv C,. The source is input on this plane by specifying complex 
amplitudes of incident duct modes [5, 10-12], The nacelle outer surface is C„. On 
this boundary, steady flow and acoustic particle velocities have a vanishing normal 
component. An artificial baffle C b formed by a ray from the origin limits the 
solution domain. The sweep angle is chosen in such a way that minimal effect on the 
acoustic field is created [13], The domain of computation is divided into two parts. 
In an inner region a standard finite-element mesh is used; in the present problem 
emht-node serendipitv elements with the condition that four to five elements per 
wavelength are required. The near field is terminated on a boundary C r beyond 
which farfield elements are used. In previous studies, this region was large but finite 
and bounded by the surface C x , a circle which represents a constant-phase surface 
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for an acoustic source located at the origin. On this boundary, a radiation 
condition was specified. Wave envelope elements [5, 10-12] were used in this 
region. In the present study, the farfield region is extended to infinity and a single 
layer of mapped infinite wave envelope elements is used to provide a reflection-free 
boundary on C r and to compute the acoustic field in the far field as required. The 
boundary C* is not part of the solution. 


3. FINITE ELEMENT FORMULATION 

The geometry of the inlet and steady flow field in and around the inlet is axially 
symmetric. The acoustic field is not axially symmetric but is represented as periodic 
in a cylindrical co-ordinate system with .v being the axis of symmetry, r the 
cylindrical radius in a circular cross-section at x = 0, and 0 the angular co-ordinate. 
Solutions are sought in angular harmonics of a Fourier series enumerated by the 
angular mode number m. This reduces the solution domain to a two-dimensional 
x, r plane. 

The starting point for the formulation of both the steady mean flow and the 
acoustic perturbation consists of the inviscid mass and momentum equations and 
the energy equation in the form of the isentropic equation of state. The acoustic 
field equations are obtained by considering small perturbations on a steady 
irrotational mean flow characterized by density p r and speed of sound c r . This 
formulation makes it possible to introduce a steady flow velocity potential (p r and 
an acoustic perturbation velocity potential (p. Acoustic perturbations in pressure, 
density and velocity potential are harmonic in time with frequency rj r and harmonic 
in the angular co-ordinate 9 of the form p(x, r)e >in, ‘~ me \ p(x, r)e‘ ( '' r ' _mB) , 
(p{x, In linearized form, the weak formulation is [5, 10-12] 




{ VW ■ (p r V(j) + p V(j) r ) - i t] r Wp) dV 


J Jv J 


f* /* 

W{p r Vcp + pV(p r )-ndS (1) 
- J 


The weighting functions are taken as W (x, r, 6) = W (x, r)e ,mfl . Angular harmonics 
proportional to e~ ime represent the decomposition of the solution periodic in 0 in 
a Fourier series. The angular mode number m is a parameter of the solution. The 
surface integral is over all surfaces bounding the domain. The unit normal for the 
surface integral is out of the domain at the surface in question. The weak 
formulation continues with the linearized momentum equation 


p = - ^(i >lr<P + V4>r- V(p) 
c; 


( 2 ) 


which is used to replace p in equation (1). The linearized equation of state. 


P = c;p. 


(3) 
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is used to produce an alternative form of the momentum equation in terms of 
acoustic pressure. 


P = - p r (ir]r<f> + V<t>r' V( P)- ( 4 ) 

The acoustic particle velocity and acoustic velocity potential are related according 
to 


v = V(p. (5) 

The linearization process also produces the weighted residual formulation for the 
steady flow, 


VW-(p r Vcp r ) dV = 


vJ 


S j 


W(p r V(j) r )-ndS, 


( 6 ) 


and the steady flow momentum equation in terms of the speed of sound. 


c? 


= 1 - 


(7 ~ 1) 


lV4>r- V<t> r - Mil 


(7) - 


and in terms of the steady flow density, 


Pr 



(7 - 1) 


1 /( 7 - 1 ) 


(V<t> r - V4>r~ Ml) 


( 8 ) 


Equations (1)— (8) are in non-dimensional form where 4> is the acoustic potential, <£ r 
is the local mean flow (reference) potential, p is the acoustic density, p r is the local 
mean flow density, p is the acoustic pressure, and c r is the local speed of sound in the 
mean flow. All quantities are made non-dimensional by using the density in the far 
field, p, o, the speed of sound in the far field, c x , and a reference length which is 
defined as the duct radius at the source plane, R, where acoustic modal amplitudes 
are defined. This plane could be the fan plane or the exit guide vane plane, but it is 
not restricted to these locations. The acoustic potential is non-dimensional with 
respect to c^R, and the acoustic pressure with respect to p x cl. Lengths are made 
non-dimensional with respect to R. Time is scaled with Rjc x , leading to the 
definition of non-dimensional frequency iy r = c oR/c x ; a> is the dimensional source 
frequency. M x = M 0 is the Mach number in the far field representing the forward 
flight effect. 

Equation (6) is the weighted residual formulation for the calculation of the 
compressible potential flow within and around the nacelle. Equations (7) and (8) 
are subsidiary relations that can be used in an iterative solution which at each stage 
uses a density field derived from the previous iteration step. V(fi r , c r , p r are required 
data for the weighted residual formulation of the acoustic problem. In the results 
reported here only the first iteration of this process is used to define the potential 
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flow field. This is accomplished by solving the incompressible problem and then 
computing a variation in steady flow density and speed of sound. 

The surface integral in equation (1) provides the boundary conditions on the duct 
walls, and at the source. The acoustic source is specified by the complex amplitudes 
of acoustic duct modes at the source plane. On this plane, the FEM nodal value of 
acoustic potential are replaced by the complex amplitudes of the acoustic potential 
modes by an eigenfunction expansion. The incident acoustic modal amplitudes are 
input and the reflected modal amplitudes are computed as part of the solution. 
Details of this procedure are available in [5, 10-12]. 

A previous study [13] shows that the baffle can be positioned to produce 
practically no effect on typical acoustic radiation patterns. Therefore, there is no 
contribution from the surface integral on the baffle. In previous studies, the surface 
integral provided the mechanism for enforcing the Sommerfeld radiation condition 
on C x . In the present application of mapped wave envelope elements the surface 
integral is never explicitly introduced on a far-field boundary because the assumed 
form of the solution in the outer region implicitly satisfies the Sommerfield 
condition. 

In terms of acoustic potential the weak formulation is, from equations ( 1) and (2), 


J v 


c ; f VW ■ V<t>-{ M r - VW){ M r - V( p) 

c: 


+ iijrlW(M r - V(fi) — (M r - VW)<p ] 


- rj;W(j)} dV 


Js 


£ [c; W V<t> - M r W ( M r • V(j)) - irj r M r W(j)} • n d S, 
c; 


( 9 ) 


where the local non-dimensional steady flow velocity is M r = V(j) r . Equation (9) is 
the weak formulation in the entire domain, however in the steady flow far field it 
simplifies considerably with the steady flow given by M r = M 0 i and p r = 1. c r = 1. 
Furthermore, the surface integral on C x has no contribution in the formulation 
proposed here because there is no longer any surface on which a Sommerfeld 
radiation condition is to be applied. The weak formulation in the steady flow far 
field is 




VW- V<p - M o tt— -r~ + irjrM 0 ( W^--^-4>] - rjfwAdV = 0. 
cx cx \ cx cx J ) 


( 10 ) 


In the cylindrical co-ordinate system used here, some liberty is taken in defining 
the gradient operations as 


cW . cW lm ccj) . c4> un , n , ,,, 

VW = i + — e r + — W^, V cp = i -t ; e r (/>e„, (11,12) 

cx cr r cx cr r 
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and suppressing factors e ±im * which arise as part of the weighting and trial 
functions as explained in connection with equation (1). These factors cancel 
throughout all of the products in equations (9) and (10). Equations (11) and (12) 
reflect the harmonic angular dependence of <j> and W. The non-dimensional 
velocity M 0 in equation (10) is the Mach number of the forward flight. 

In the steady flow near field, where the flow is non-uniform, equation (9) is 
discretized using standard finite-element techniques. Example calculations 
presented in this study are based on two-dimensional rectangular isoparametric 
serendipity elements with eight nodes. 

In the steady flow far field where the flow is essentially uniform, equation (10) can 
be discretized using wave envelope elements or by. introducing mapped infinite 
wave envelope elements to obtain closure of the computational domain. It is the 
formulation in terms of mapped wave envelope elements which is of interest here. 


4. THE INFINITE MAPPING 

Because of the harmonic dependence on the angle 0 the originally 
three-dimensional weak formulation reduces to two spatial co-ordinates x and r. 
The x, r plane is shown in Figure 1 where the boundary C r separates an interior 
region in which standard FEM descretization is used from an outer region in which 
mapped wave envelope elements are used. The exterior region must be in the steady 
flow far field. Figure 2 shows an element in the outer region in the x-r plane of the 
cylindrical co-ordinate system. The element is bounded by the surface C r on which 



x 


Figure 2. Details of the finite/infinite-element interface. 
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Rays 



Figure 3 . Geometric details of the mapping between the infinite element and the parent element: (a) 
parent element, (b) mapped wave envelope (WE) element. 


it conform with an element of the conventional mesh. The edges of the element are 
straight lines, extending outward more or less radially, though not necessarily from 
the axis system origin nor necessarily from a common origin. For the elements used 
in this investigation which conform with eight-node serendipity elements in the 
conventional mesh (each with three nodes on C r ), a third radial line between the 
two edges is required. For simplicity, each of the three straight lines will be referred 
to as rays. In Figure 1, the outer surface C x is the notional outer boundary of the 
element at infinity. A ray of an element has an apparent origin at a point x 0 ,r 0 
which in general can be different for each ray. The element maps to a parent 
element in the c, /? plane. — 1 < q ^ 1, — 1 < n < 1, as shown in Figure 3. The rays 
of the element map to the c axis with q = — 1, 0, 1 in the parent element according 
to 


Since 





1 + c 



(13, 14) 



1 C 


i. 


(15) 


the mapping is unchanged by an origin shift. Therefore, it can also be used to yield 
a mapping relative to the source at x 0 , r 0 : 


x - x 0 



— - x o) + 3 — ' — ~(- x 2 — - x o)> 


(16) 
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r ~ r o = y—:(ri ~ r 0 ) + {-^(r, - r 0 ). (17) 

The node x L , r t is defined by the requirement that the element conform with the 
conventional element on the boundary C r . A particularly useful form of the 
mapping is obtained if the node at x 2 , r 2 is located such that x 2 — x 0 = 2(.v , — x 0 ) 
and r 2 - r 0 = 2(r x — r 0 ). This makes the mapping simplify to 


X - Xn = 


2( -Vl - -To) 

1 - c ’ 


/*0 — 


2( r L -r 0 ) 

1 - C ‘ 


(18) 


The mapping has the properties that <; = — 1 maps to .v — ,v 0 = xq — x 0 . 
r — r 0 = r t — r 0 , c = 0 maps to x — x 0 = 2 (.t, — x 0 ), r - r 0 = 2 (r x — r 0 ), and 
c = 1 maps to x — x 0 = x , r — r Q = x. The mapping along a “ray" transforms the 
infinite domain in the physical co-ordinates to the domain — 1 ^ c 1 in the 
parent element. The inverse mapping is 


2Ul - .to) 



2(r t - r 0 ) 
r ~ r 0 


(19) 


It is easily deduced that this mapping along a “ray" also applies for the polar 
radius of the point x, r relative to x 0 , r 0 , r p = v /(x — x 0 ) 2 + (r — r 0 ) : , in the 
form 


2 v '(.t 1 - -T 0 )~ + (r, - r 0 ) 2 

r„ =— : : 



( 20 ) 


and 

i-l-2rjr r . CD 

This form emphasizes the role of the base node x 0 , r 0 as a “source” for the “ray" and 
the distance r p as the polar distance from the source. 

The infinite-element mapping is completed by a conventional mapping on 
— 1 ^ rj ^ 1. The element shown in Figure 3 has six nodes numbered as shown. 
Nodes 1-4 are corner nodes and nodes 5 and 6 are mid-side nodes on C r and C 2 
(the locus of the nodes x 2 , r 2 ). The mapping is of the form 


x = [M(c, > 7 )]x, r = [M(c ? /?)]i\ 


where [M (c, ^)] is a row vector of six shape functions AT(T rj) and x, r are vectors of 
nodal values of x, r. With the nodal numbering scheme used in Figure 3 the explicit 
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M { (c, n) = 0-5)l(>l - 1)- — — : 

L — C 


= 05 



_ 1 - 

rj) = 0-5 /;(// 4- 1)- — ^ 

i — L 


iV/ 4 (s, n) = 05/7 ('/ 



M 5 (C,/7) = (1 -//)(1 +/7)-^. 

i L. 


M b (c. //) = (1 - /?)(! + n) 


i + ; 
1 - c 


The mapping described here is simply another view ot the results presented by 
Astley et al. [8] specialized to the cylindrical co-ordinate system. 

In preparation for development of a mapped infinite wave envelope element for 
a uniformly flowing medium it can also be noted that the results of 
equations ( 1 6)— (2 1 ) can be extended to other "distances" along a ray yielding 
a similar mapping. For example 

2R 

R = v /(- x ~ - x o) 2 + P 2 (r ~ r 0 ) 2 = 7 4- <- 4 ) 

1 — c 


where R i = ^/(x, - x 0 ) 2 + /F(r, - r a ) 2 and p 2 = l - M 2 and M is the Mach 
number of the uniformly flowing medium. A second useful mapping is 

= i[- M(x - x 0 ) + R]= (25) 

P~ 1 ~ 

where ^ = (1//J 2 )[ — A/(Xi - x 0 ) + K t ]. These observations are important to the 
extension of the application of the mapped infinite wave envelope element to 
acoustic radiation in uniform steady flow. 


5. SOURCE SOLUTION IN UNIFORM FLOW 

The weak formulation of equation (10) for acoustic radiation in a uniformly 
moving medium is consistent with the differential equation 


7T + jV/tt-V </> = V 2 (j). 
ct cx } 


(26) 


A fundamental harmonic source solution for this equation is 


0 


a i nrt 


e ( - - V/.x - v .X- + li-r- ) 




(27) 
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where /? = V T — M~ and r 2 = y 2 + r 2 . This can be verified by direct substitution 
or by noting that the transformation of variables 

*-l K-r. + «» 

reduces the convected wave equation (26) to the standard wave equation 

5 2 <i>/dt’ 2 =V 2 cj) (29) 

in the transformed variables which has a fundamental harmonic source solution 


4> 


.Hi . e -i /^ * "~ r 

= e 1 "' ,, , 


+ r 


(30) 


Equation (27) is then obtained replacing the change of variables of equation (28). In 
terms of the definitions of equations (24) and (25), the fundamental harmonic 
source solution for source location at ,x 0 , r 0 is 

p. - 

</> = e ir, '‘ — - — . (31) 

R 


6 . SHAPE FUNCTIONS IN THE INFINITE ELEMENTS 

Shape function in the mapped infinite wave envelope elements can be 
constructed to display the characteristics of the fundamental source solution at 
large distances from the source in the form 

p ~ i rj,OMx) -»//,) 

<t> = CMx)e~ im<) /?, — = Pfxle-^e-^", (32) 

K (x) 

where the notation x = (.x, r) and j.i (x) = ry r (i^(x) — i/'i) is used and is similar to the 
notation used by Astley et al. [8]. /t(x) is the phase relative to the surface C r 
separating the infinite-element region from the region of standard FEM 
interpolation and i/q emphasizes that this phase is dependent on the specific ray 
on which equation (32) is evaluated, i/q would be a constant for the entire 
infinite-element region if C r is a surface of constant \p (a constant-phase surface ), 
but in general would vary from node to node on C r . The most direct way to make 
i/q invariant for the mesh is to construct the mesh so that for all infinite elements 
-To< r o (the “source point”) is common and C r is a surface of constant tjy ( phase ) 
relative to the common “source”. The mesh used by Eversman et al. [5, 10-12] has 
this property (.x 0 , r 0 are at the mesh origin) and is used in examples in this 
investigation. At large R, equation (32) should have asymptotic behavior in 
R consistent with equation (31). The function P(x) should therefore display the 
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Rays 



Figure 4. Example of an infinite element with nine x interpolation nodes and six • mapping nodes. 
This element produces an asymptotic interpolation in the far field of third order in Ri/R. 


appropriate asymptotic behavior in R, should be capable of accounting for 
nearfield effects, and should interpolate in the standard FEM context in the 
rj co-ordinate in the parent element. 

In terms of <j, rj co-ordinates of the parent element, ju(x) and x , have simple 

forms suggested by equations (24) and (25); 




r(c, n) = <Ai 


i±i Ri 

l-c’ R(c,rj) 



C). 


(34, 35) 


In equation (34), i/q can be a function of 17 on the inner boundary of the element 
c = — 1 , interpolated relative to nodal values on C r . The function 0(x) in 
equation (32) which accounts for nearfield behavior in the infinite element can be 
represented by a standard FEM interpolation 

Q(l n) = [S(c, rj)10_, (36) 

where Q is a vector of nodal values of (?(x). There are six nodes involved in the 
infinite mapping and these can be used as nodes in the interpolation of Q(x). It will 
generally be appropriate to use more than the mapping nodes by including extra 
nodes along the rays as shown in Figure 4 which demonstrates the introduction of 
one extra node midway between the mapping nodes on each ray and suggests 
a convenient nodal numbering scheme. 

The shape functions for the element shown in Figure 4 with the additional node 
midway between the mapping nodes on a given ray are based on nine-node 
Lagrangian interpolation with the extra nodes mapped to c = — 3. In general, P(x) 


38 



MAPPED INFINITE ELEMENTS 


677 


is interpolated within an element according to 


P{c,r 1 ) = lN(Q,r 1 )}Q, (37) 

where the shape function N,{q, rj), the shape function corresponding to node i, is 
constructed from the mth-order Lagrangian shape function for node i, Lfd, rj), 
according to 


p,(c, n) = NAi, n) = i(i - n )■ (38) 

Some liberty is taken with notation here; Lf(c, 17 ) is defined so that p is the order 
of interpolation (number of nodes) along the c-axis. Along the / 7 -axis the order 
conforms with the order used in the standard FEM region, which is 3 in the 
two-dimensional serendipity elements implemented in the model reported here. It is 
interesting to note that N)(<;, rj) is unity only for nodes with c, — — 1 (the Lagran- 
gian interpolation functions have the value unity for all nodes). The vector of nodal 
values Q corresponding to the evaluation of <2(x) only corresponds to nodal value 
of P(x) for nodes on the surface C r . Because of this, and because of the phase term 
e - '" (x) in equation (32), which is unity only on the surface C r , in the infinite elements 
the solution vector does not correspond to nodal values of acoustic potential at 
most of the nodes. The potential can be easily reconstructed by postprocessing. 

The form of the shape functions defined by equation (38) can be interpreted in 
global co-ordinate by using equation (24) to show that 

1 - { = 2 (RJR), d = 1 - 2 (RJR). (39) 

Equations (38) and (39) suggest that the shape functions in global co-ordinates 
along a ray are of the form 

P,(x, r) = y t (RJR) + y 2 (Ri/R) 2 + 7i(Ri/R) 3 + - + y.(*iW. ( 40 > 

n is determined from the order of Lagrangian interpolation. For a p node 
interpolation leading to polynomials in q of degree p — 1 it is determined that n = p. 
A similar result was shown in the case of radiation in a stationary medium [ 8 ]. 

Reference to “variable order” mapped infinite wave envelope elements relates to 
the choice of the order of the Lagrangian interpolation and therefore to the powers 
of RJR in the asymptotic expansion for the shape function. Conceptually this 
could be extended to any order, but as pointed out by Astley et al. [ 8 ] there is 
a limit imposed by the onset of numerical problems probably related to ill 
conditioning if the order is too high. 

7. WEIGHT FUNCTIONS IN THE INFINITE ELEMENTS 

Astley et al. [ 8 ] show that in order for the boundary integral introduced in the 
weak formulation to have no contribution on the boundary at infinity it is 
necessary for the weighting functions to be functions of {Ri/R{x)} q+1 , with q > 1. 
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The weight functions are of the form 


</> = 0(x)e imd 



_ / )(\)P(x)e imtf e^ (xl 


(41) 


where 


D(x) = (Ri/RWy 


(42) 


In the parent element, 


D(c. n) = (lz-Hl -cY (43) 

The weight functions are the complex conjugates of the shape functions multi- 
plied by the additional decay term. In the present investigation q = 3. The notation 
here has been chosen to correspond to that used by Astley et al. [8] to emphasize 
the similarity with their development in the case of a stationary medium. Only the 
details hidden in the definitions of i p and R are different. 


8. THE WEAK FORMULATION IN THE INFINITE-ELEMENT REGION 

The weak formulation of equation (10) for the infinite-element region in which 
the steady flow is necessarily uniform is obtained by using equations (32) and (43) 
defining the assumed form of solution and the weight functions in the infinite- 
element region. The gradient operations on the assumed shape and weighting 
functions yield 


and 


V<f> = ( VP - iiiT^fe"" 1 '" 


(44) 


VW = (DVP* + ii/ r DPF )< + PVD)e Uh>1 . 


(45) 


where the notation 


VP 


cP . cP i m 

— i + — e r Pe 0 , 

cx cr r 


VP* 


cP . cP i in _ 

— — 1 H e r 4 

cx cr r 


(46. 47) 


is used as in equations (1 1) and (121) because of the factors e ± ' m0 which are 
suppressed. By using standard finite-element operations, equation (11) can^be 
formulated at the global level to yield complex element "stiffness” matrices [K, ; ] 
defined in terms of real mass, stiffness and damping matrices, 

[K„] = - viMul + iG[C, ; ] + [K„], (48) 


40 


MAPPED tNFINITE ELEMENTS 


679 


where 
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(49) 

(50) 


C U = 
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J J 

V j 

1 
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,, CLIO Pi ,,cPj CficPi 

D (1 - + M-z 1 + C~ — 

cx cx cx cr cr 


cD 3jj. ,3D dD CfjL 

(1 -M 1 ) — Jt + M— + — -f 
cx cx ox cr cr 


PiPj> dV. 


(51) 


The definitions of the stiffness, mass, and damping matrices of 
equations (49)— (5 1) are implemented at the element level using the infinite mapping 
to the parent element. These results reduce to those of Astley et al. [8] when the 
medium is stationary and when account is taken of the operations which are 
particular to the cylindrical co-ordinate system. It is not difficult to generalize to 
a three-dimensional Cartesian co-ordinate system. 


9. AN IMPORTANT PROPERTY OF THE MASS MATRIX 

The mass matrix of equation (50) vanishes if the surface C r separating the 
standard finite-element region from the infinite-element region is a surface of 
constant phase for an apparent acoustic source location Xq, r q which is common for 
all elements. This is shown by referring to the definition of ju(x), 


fl(\) = *lr( |A(X) - 


(52) 


where 


i/Tx) = (l/j 3 2 )[- M(x - x 0 ) + K], R = y/(x - x 0 ) 2 + j3 2 (r -r 0 ) 2 . (53,54) 

Since it is stipulated that C r is a constant phase surface, it follows that i/q is 
constant. The apparent source location is the same for all elements, leading to the 
conclusion that x 0 , r 0 are constants. It can then be verified that 


(1 - M 2 )(du/8x) 2 + [on! dr) 2 + 2 M(dfi/dx) = 1 


(55) 



680 W. EVERSMAN 

which from equation (50) leads to the result 

= 0. (56) 

This is consistent with the findings of Astley et al. [8] in the case of a stationary 
medium when the surface C r is a sphere, a constant-phase surface in this case. While 
of some interest in the time harmonic formulation considered here, the vanishing of 
the mass matrix is of central importance when a time-dependent formulation is 
implemented in the stationary medium case. It remains to be established that this is 
equally important in the case of a uniformly moving medium. 

10. TURBOFAN INLET EXAMPLE 

Figure 1 shows the generic geometry of a turbofan inlet in an x, r plane of 
a cylindrical co-ordinate system. The nacelle interior and exterior shape are typical 
of realistic nacelles. The acoustic source is on the plane C ' j and produces a combi- 
nation of radial modes at a specified angular mode m and non-dimensional 
frequency rj r . The source strength is specified by the complex mode amplitudes. 
This type of source would correspond to rotor alone noise or rotor/exit guide vane 
interaction noise. The frequency is determined by the number of blades on the rotor 
and the angular mode number by the rotor and exit guide vane blade counts. The 
nacelle has a forward velocity specified by the Mach number M 0 , which is 
represented for the stationary nacelle by a steady flow directed toward the nacelle. 
The steady flow into the nacelle is specified by the Mach number taken to be 
uniform on the source plane. The steady flow field inside and outside the nacelle, 
computed on the FEM acoustic mesh, provides input data for the FEM acoustic 
calculations. This mesh is over refined for the steady flow calculations but this 
inefficiency is more than offset by the convenience of input data on a mesh 
compatible with the acoustic mesh. 

The details of the FEM acoustic computations with the domain closed by 
a conventional wave envelope transition region to a Sommerfeld radiation bound- 
ary are given in references [5, 10-12]. In this example the propagation and 
radiation problem is formulated with the standard FEM treatment in the steady 
flow near field and the domain is closed in the far field by the use of mapped infinite 
wave envelope elements. The specific case shown is at a reduced frequency rj r = 25 
and angular mode m — 23 with only the first radial mode incident. Only one radial 
mode propagates and it has a cutoff ratio near unity, which indicates that the 
peak lobe of the radiation pattern will be at a high angle relative to the nacelle axis. 
In this case it is over 60° to the nacelle axis for the case of M, = 0-20 and Af 0 = 0-30. 
Figure 5 shows the standard mesh in the region which has been abritrarily declared 
as the steady flow near field. The steady flow far field is where the flow is essentially 
the jV/ 0 = 0-3 uniform flow. The outer boundary of this mesh is the surface C r and it 
is a circle of constant phase for a source at the axis system origin. The infinite- 
element region is outside of C r and not shown. The same inner mesh was used 
with the outer region consisting of seven layers of standard wave envelope 
elements extending to 10 duct radii ahead of the inlet for the purpose of producing 
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Figure 5. The nearfieid mesh of standard finite elements bounded by the surface C r . 


15.0 T 



Fisure 6. Standard wave envelope mesh used in finite element/wave envelope element formulation. 

comparison results. The outer mesh for this case is shown in Figure 6. The standard 
code has been extensively benchmarked by experiment [12] and by comparison 
with available approximate analytical results. Numerical experiments have 
shown that for this frequency radiated fields are particularly difficult to model. 
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Figure 7. Contours of equal acoustic potential in the entire computational domain for the finite- 
element/wave envelope element formulation. External mach number M 0 = 03, source plane Mach 
number A/, = 0-2, non-dimensional frequency rj r = 25, angular mode number m — 23, first radial 
mode. 


At high angles to the axis the source is certainly not seen as a simple source on C r as 
located in this example. It is reasonable to expect that non-reflecting boundary 
behavior based on an asymptotic approximation representing a simple source 
would be difficult to achieve. 

The results which will be displayed are contours of constant acoustic potential 
magnitude in an .x, r plane superposed on the nacelle geometry. Acoustic potential 
has been chosen since there is an extra post-processing step to obtain acoustic 
pressure which introduces its own potential for error, unrelated to the details of the 
reflection free boundary. Post-processing for pressure in the standard FEM region 
involve the same operations whether standard wave envelope or mapped infinite 
wave envelope elements are used in the outer solution. Figure 7 shows the radiation 
pattern generated by using the standard code (wave envelope elements) and 
plotting contours of constant acoustic potential in the entire computational do- 
main. Figure 8 shows the same results limited to the region of standard finite 
elements, which provides a more detailed way of viewing the reflection free perfor- 
mance of the boundary C r . Figure 9 shows the results when mapped infinite wave 
envelope elements are used to provide a reflection free boundary. In the case shown 
the formulation is based on eight-node Lagrangian interpolation in the mapped 
elements in the c direction (eight nodes). This corresponds to introducing R\/R in 
the expansion for asymptotic behavior of the farfield solution up to the eighth 
power [refer to equation (40)]. Element integration is based on 9 x 3 Gauss points. 
It was found that five-node Lagrangian interpolation (powers of R[/R up to five in 
the asymptotic expansion) was not sufficient. 
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Figure 8. Contours of equal acoustic potential in the standard finite-element region for the 
finite-element/wave envelope element formulation. External mach number iV / 0 = 0'3, source plane 
Mach number j V/, = 0-2. non-dimensional frequency rj r = 25, angular mode number m = 23, first 
radial mode. 


Figure 7 displaying the entire solution field to the Sommerfeld boundary (10 
ducFradii of the nacelle axis) suggests significant reflection from the boundary 
which appears in the waviness of the contours, particularly at higher angles where 
diffraction around the inlet lip is important and where the nacelle surface interferes 
with the radiation. The quality of the solution does not improve with the further 
mesh refinement, indicating that the mesh is suitable for the frequency. Figure 8 
zooms in on the region inside C r and the poor quality of the solution is apparent. In 
Figure 9 the same level contours are considerably less ragged, indicating that 
reflection has been essentially eliminated. It is of interest to recall that the computa- 
tional domain includes the artificial baffle C b and it appears that it has little effect 
on the radiated field, consistent with the results reported in reference [13]. 

The clear conclusion is that poor quality of the solution when standard wave 
envelope elements are used is due to the inability of the wave envelope elements to 
provide a completely reflection-free boundary for the complicated source config- 
uration and this location of C r . In principle expanding C r should improve the wave 
envelope element performance, but this has the obvious implication ot directly 
increasing the dimensionality (presuming it is required that the mesh refinement is 
retained) and the hidden implication of requiring even further mesh refinement due 
to the growth in element aspect ratio as C r is expanded. 

Variable order mapped infinite wave envelope elements generally will increase 
the maximum bandwidth of the mesh (the inner mesh may have eight nodes per 
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Figure 9. Contours of equal acoustic potential in the standard finite-element region for the 
finite-element/mapped infinite wave envelope element formulation. External mach number M 0 = 0-3. 
source plane Mach number A/, = 0-2, non-dimensional frequency rj r = 25, angular mode number 
m = 23, first radial mode. The order of the interpolation in the infinite elements is eighth power in 
RJR. 


element and the single infinite element layer has been tested here with as many as 24 
nodes per element). In frontal solvers, this tends to slow down the solution even if 
the total number of nodes is more or less the same in the mapped elements as in the 
standard wave envelope elements. In the example discussed, a frontal solver is used 
and the mapped infinite-element computation have an execution time which is in 
a ratio of about 7/5 compared to the standard wave envelope code. This cost is not 
unimportant, but must be assessed against the requirements for solution quality. In 
this case, the infinite-element results are clearly superior. 

The question now arises; how much can the computational domain be reduced 
by using the infinite elements to enhance the reflection-free boundary? To partly 
address the question, the boundary C r has been reduced to a radius of two duct 
radii ahead of the origin. Note in the original mesh of Figure 5 the mesh extends 2-5 
duct radii ahead of the origin. In order to maintain approximately the same 
mesh refinement, the element count between the “highlight circle - ’ (a circle 
passing through the tip of the inlet lip and intersecting the axis near r = 1) has 
been reduced from 50 to 35. Figure 10 shows acoustic potential level curves in 
the standard element region for the case using mapped wave envelope elements 
for closure. The quality of the solution is still substantially superior to that seen in 
Figure 8 for which closure was achieved using regular wave envelope elements 
(note that the level curves are not the same in Figures 8 and 10 because they are 
based on the maximum level on C r , which differs because C r differs). The 
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X/A 

Figure 10. Contours of equal acoustic potential in the standard finite-element region for the 
finite-element/mapped infinite wave envelope element formulation on a reduced mesh in the standard 
element region. External mach number M 0 = 0-3, source plane Mach number M, = 0-2, non-dimen- 
sional frequency rj r = 25, angular mode number m = 23, first radial mode. The order of the interpola- 
tion in the infinite elements is eighth power in Ri/R. 


computation time ratio is now nearly 1/1 and the mapped infinite-element results 
are still superior. 


11. CONCLUSION 

It has been shown that with suitable modifications mapped infinite wave envel- 
ope elements can be used to provide an effective reflection-free boundary for 
acoustic radiation in a uniform steady flow. The adaptation of the elements to this 
case is based on the observation that all important “distances” along “rays” map to 
the parent element in the infinite mapping in exactly the same way. This permits the 
fundamental solution for radiation from a source in uniform flow to be mapped to 
the parent element in a form similar to the mapping in the case of a stationary 
medium. The fundamental solution forms the basis for an asymptotic expansion in 
R~ n in the infinite elements, where R is the “convected radius”, R 2 = x 2 + fi 2 r 2 . 
The order of the asymptotic expansion can be chosen to meet the needs of the 
problem. Element mapping functions are identical to those previously proposed for 
the stationary medium case and the shape functions are of the same form as those in 
the stationary medium case with differences only in the details. 

Computational examples have been based on acoustic radiation from a turbofan 
inlet which has been the subject of several previous investigations in which an FEM 
model was developed with the reflection-free closure of the computational domain 
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based on standard wave envelope elements. Examples have shown that mapped 
infinite wave envelope elements provide a superior reflection-free boundary for 
cases in which the standard wave envelope elements generate reflections which 
appear in the radiated field. It should be noted that the improved performance may 
not be without cost. If relatively high order mapped elements (asymptotic behavior 
to R~ n where n is relatively large) are required, the maximum front width of the 
FEM formulation may be larger than would occur in the standard wave envelope 
element formulation. For frontal solvers this may decrease computational efficien- 
cy. However, this cost has a substantial benefit in the quality of the solution which 
may not be achievable with the standard wave envelope elemnts without expanding 
the boundary between standard FEM and the wave envelope element region. In 
fact, it has been shown that by taking advantage of the reduction in size of the inner 
region (standard element region) which is achievable with mapped infinite elements 
it is possible to obtain superior solutions without increasing computation time. 

It has been found that the mapped infinite wave envelope element region can be 
cast in the form of appended mass, damping and stiffness matrices. With a suitable 
choice of the surface which separates the standard FEM region from the infinite- 
element region and the restriction that the mapping and shape functions in the 
infinite elements are based on a common apparent source location, it has been 
shown that the element mass matrices vanishe. This has previously been shown to 
be important for transient FEM formulations for radiation in a stationary medium. 
This suggests that similar investigations should be carried out in the case of 
uniform external flow. 
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ABSTRACT 

Variable order mapped infinite wave envelope elements are developed for finite element 
modeling of acoustic radiation in a uniformly moving medium. These elements are used as a non- 
reflecting boundary condition for computations on an infinite domain in which a radiating body 
is immersed in a moving medium which is essentially undisturbed outside of the near field. The 
mapped elements provide a boundary condition equivalent to element stiffness, mass, and damping 
matrices appended to an inner standard FEM mesh. A demonstration of the performance of 
mapped elements as influenced by element order is given in the context of acoustic radiation from 
a turbofan inlet and exhaust. 

INTRODUCTION 

FEM modeling of acoustic radiation is usually complicated by the requirement that 
prediction of the acoustic field is sought in some finite sub-domain of an infinite domain. This 
dictates that computations be limited to the sub-domain with a non-reflecting boundary or that 
the infinite domain be mapped to a finite computational domain. The investigation reported here 
seeks an improved non-reflecting boundary condition for aeroacoustics problems in which 
acoustic radiation in the far field is influenced by steady uniform flow. In order to maintain 
compatibility with standard FEM solution procedures (for example frontal solvers or other sparse 
matrix/narrow bandwidth methods), only local types of boundary conditions are considered. In 
this context boundary conditions which can be categorized as various evolutions of “infinite 
element” methods [1-5] or “wave envelope element” methods [6,7] have been studied. Mapped 
infinite wave envelope elements [4,5], which combine the attributes of “infinite” and “wave 
envelope” elements, have been chosen in this investigation because of their almost seamless 
compatibility with active FEM codes and meshes. Mapped infinite wave envelope elements limit 
computations to a finite domain and provide an approximate reflection free boundary. They have 
been investigated extensively for acoustic radiation in a stationary medium [4,5], Formulation of 
the elements reveals the possibility of including within the element shape function an explicit, and 
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adjustable, dependence on inverse powers of the distance from an apparent acoustic source. This 
allows the introduction of mapped infinite wave envelope elements well into what would normally 
be considered the acoustic near field, reducing mesh refinement and dimensionality. 

The study reported here extends the variable order mapped infinite wave envelope element 
concept, as implemented by Astley and co-workers [4,5], to steady uniform flows, principally in 
connection with aeroacoustic problems related to turbofan acoustic radiation. Most previously 
reported applications of infinite elements were implemented in two or three dimensions. The 
application here is in a cylindrical coordinate system reduced to two dimensions by taking 
advantage of periodicity of the solution in the angular coordinate. References [3-5], and a recent 
contribution by Astley and Hamilton [8], provide excellent citations to the development of the 
infinite element concept, leading to the present application. 


AN APPLICATION TO TURBOFAN ACOUSTICS 

An important problem of acoustic radiation in a moving medium is available in the study 
of the acoustic field of a turbofan inlet and exhaust. Acoustic propagation and radiation occurs 
in a high speed potential flow which is the net effect of flow into (out of) the inlet (exhaust) and 
the forward flight of the nacelle. In the steady flow far field acoustic radiation occurs in a 
unif ormly moving medium. It is required to make computations to predict the radiated field 
relatively near the inlet [6] (exhaust [7]). This has been approached in the past by terminating the 
computational domain with a Sommerfeld condition on a boundary reached by the use of wave 
envelope elements. Here it is intended to investigate the application of mapped infinite wave 
envelope elements to obtain closure of the computational domain. 

The problem is cast here in terms of the turbofan inlet, but an exhaust flow can be modeled 
with modifications particular to the shear layer boundary between the exhaust jet and the 
surrounding medium [7], The geometry of the nacelle and steady flow field in and around it is 
axially symmetric. The acoustic field is not axially symmetric but is represented as periodic in a 
cylindrical coordinate system with x being the axis of symmetry, r the cylindrical radius in a 
circular cross section at x = 0 , and 0 the angular coordinate. Solutions are sought in angular 
harmonics of a Fourier Series enumerated by the angular mode number m . This reduces the 
solution domain to a two dimensional x , r plane, shown in Figure 1. The inlet shape in a 
0 = constant plane is defined by the surface C n which includes the center body. The surface C- 
is the plane on which a source is defined, for example the plane of the fan. The surface C b is an 
artificial baffle introduced to limit the computational domain. The boundary C„ is the outer 
boundary of the computational domain, which in principle is infinitely far away, but may be a finite 
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surface where a radiation condition is introduced. 

The acoustic field is assumed to be harmonic in time at non-dimensional frequency r) r 
Geometry is non-dimensional based on a reference length generally chosen as the radius of the 
inlet at the source plane, R . Acoustic and steady flow variables are non-dimensional based on 
reference values of the speed of sound and density of the medium, p_ , c., generally defined in 
the uniform exterior flow. The non-dimensional frequency is q r = a) R / c m , with cnthe harmonic 
source frequency. 

In terms of acoustic potential the weak formulation is [3] 



where the local non-dimensional steady flow velocity is M r = Vcj) r and the local non-dimensional 
density and speed of sound are p r , c r . The surface integral on the right hand side introduces the 
noise source on C f in Figure 1 and a possible impedance boundary condition on C n inside the 
inlet. Equation (1) is the weak formulation in the entire domain, however in the steady flow far 
field it simplifies considerably with the steady flow given by M r = M 0 i and p r = 1 , c r = 1 . 
The weak formulation in the steady flow far field where the flow field is uniform at Mach number M 0 
is 


f f f {VW’Vfy - Ml — ^ -i\ 2 r W$}dV 

J J J dx 3x 3x OX 

V 

= II W{V<b - + ir\M^)T}-ndS 

S 


( 2 ) 


Equation (2) is obtained from equation (1) by noting that M r is constant and directed along the 
x axis. The surface integral on the right is only on the outer boundary C„ and provides the 
possibility of introducing the Sommerfeld condition, if required. 

In the steady flow near field, where the flow is non-uniform, equation (1) is discretized 
using standard finite element techniques. Example calculations presented in this study are based 
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on two-dimensional rectangular isoparametric serendipity elements with eight nodes. In the steady 
flow far field where the flow is essentially uniform, equation (2) is discretized by introducing 
mapped i nfin ite wave envelope elements to obtain closure of the computational domain. 

THE INFINITE MAPPING 

Only a brief summary of the infinite mapping is given here, as it is given in detail elsewhere 
[4], Equation (1) is discretized using standard FEM techniques in a near field region bounded by 
the curve C r shown in Figure 1. In a far field region, bounded by C r and C.,a notional 
boundary at infin ity equation (2) is discretized using mapped infinite wave envelope elements. 
These elements conform with standard elements in the inner standard FEM region on C r . Figure 
2 shows an infinite element in the outer region. The edges of the element are straight lines, 
extending outward more or less radially, though not necessarily from the axis system origin nor 
necessarily from a common origin. For the elements used in this investigation a third radial line 
between the two edges is required. A “ray” of an element has an apparent origin at a point 
x o > r o which in general can be different for each ray. The element maps to a parent element in the 
E, , T) plane, -l<l;<l,-l<q<l as shown in Figure 2. The rays of the element map 'to 
the E, axis with q = - 1 , 0 , 1 in the parent element according to 


x ~ x o 


2{x l - x 0 ) 
1 -§ 


r - r o = 


2(r t - r 0 ) 

1 -S 


( 3 ) 


It is easily deduced that this mapping along a “ray” also applies for the polar radius of the point 
x , r relative to x 0 , r Q , r p = ^{x - q,) 2 + (r - r 0 ) 2 , in the same form 


2 y / (x 1 - x 0 ) 2 + (r t - r 0 ) 2 _ (4) 

r? = 1 - § " 1 - 5 

This form emphasizes the role of the base node x 0 , r 0 as a “source” for the “ray” and the 
distance r as the polar distance from the source. The infinite element mapping is completed by 
a conventional mapping on -1 < q < 1 . It is also noted that the results of equation (4) can be 
extended to other “distances 11 along a ray yielding a similar mapping. For example 

R = f x " * 0 ) 2 + PV - ''o ) 2 = J-J ^ 

where R x = ^ - x a ) 2 7 and p 2 = 1 - M 2 and M is the Mach number 
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of the uniformly flowing medium. A second useful mapping is 

i|f = - r 0 ) + /?] = (6) 

where iJj, = P' 2 [-M(x t - x 0 ) + /?J. These observations are important to the extension of the 

application of the mapped infinite wave envelope element to acoustic radiation in uniform steady 
flow. 

The weak formulation of equation (2) for acoustic radiation in a uniformly moving medium 
is consistent with the differential equation 


( — + A/— ) 2 4> = V 2 c}) 
3t 3x 


(7) 


Equation (7) has a fundamental harmonic source solution at non-dimensional frequency r| r , for 
source location at x 0 , r 0 , given by 


c£> = e 


<V e 


'Vf 1 


R 


(3) 


where i|r and R are defined by equations (5) and (6). Equation (8) can be deduced by 
transformation of equation (7) by noting that the transformation of variables 


x' = £ , r' = r , t' = + pf 


P ' ' P 

reduces the convected wave equation to the standard wave equation 


(9) 


C 2 (J) 


= v /2 4> 


dt 


/2 


( 10 ) 


Equation (10) in the transformed variables has a fundamental harmonic source solution at 
frequency r| / r = r\ r / P which is 


{7 r Tr n 

Equation (8) is then obtained (within a constant) by replacing the change of variables of equation 
(9) and accounting for the source location at x Q , y Q , z 0 . 
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SHAPE FUNCTIONS IN THE INFINITE ELEMENTS 

Shape functions in the mapped elements can be constructed to display the characteristics 
of the fundamental source solution at large distances from the source in the form 


“'TlrW*) “ W 

<b = 0(x)e ' /mS R. = P(x) e ~"" 0 e 0 2 ) 

y ~ 1 R(x) 

where the notation x = (x , r) and p(x) = q r (t|r(x) - ^i) is used. p(x) is the phase relative 
to the surface C r separating the infinite element region from the standard FEM and ijr, 
emphasizes that this phase is dependent on the specific “ray” on which equation (12) is evaluated. 
tJ/j would be a constant for the entire infinite element region if C r is a surface of constant ij/ (a 
“constant phase surface”). The most direct way to make ijq invariant for the mesh is to construct 
the mesh so that for all infinite elements x Q , r 0 (the “source point”) is common and C r is a 
surface of constant r|r (“phase”) relative to the common source . The mesh used by Danda Roy 
and Eversman [6] has this property (x 0 , r 0 are at the mesh origin) and is used in examples in 
this investigation. At large R equation (12) should have asymptotic behavior in R consistent 
with equation (8). The function P(x) should therefore display the appropriate asymptotic 
behavior in R , should be capable of accounting for near field effects, and should interpolate in 
the standard FEM context in the q coordinate in the parent element. 

In terms of the E, , q coordinates of the parent element p(x) and R(x)/ R { have 
simple forms suggested by equations (5) and (6): 

g(5 , H) * ^ ( 13 ) 


gi 

R(Z, , q) 


= ^(1 -0 
2 


(14) 


In equation (12) in the most general case ijq can be a function of q on the inner boundary of 
the element £ = - 1 , interpolated relative to nodal values on C r . 

The shape function P t (^ , q) corresponding to node i is constructed from the p th order 
Lagrangian shape function L? (Z, , q) , according to 
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( 15 ) 


P t a, ti) =^(^,7!) = 1(1 -5)V(5.n) 

Some liberty is taken with notation here; Z'(5 , n)« defined so that p is the order of 
interpolation (number of nodes) along the 5 axis. Along the rj axis the order conforms with the 
order used in the standard FEM region, which is 3 in the two dimensional serendipity elements 
implemented in the model reported here. 

The form of the shape functions defined by equation (15) can be interpreted in global 
coordinates by using equation (1 5) to suggest that the shape functions in global coordinates along 
a ray are of the form 


R, R u . 
Y!(Y> + Y 2 (y ) 2 




(16) 


n is determined from the order of Lagrangian interpolation. For a p node interpolation leading 
to polynomials in \ of degree p - 1 it is determined that n - p . A similar result was shown 
in the case of radiation in a stationary medium [4], The mapped elements are variable order 
because n can be chosen for the application. 


WEIGHT FUNCTIONS IN THE INFINITE ELEMENTS 

In order for the boundary integral introduced in the weak formulation to have no 
contribution on the boundary at infinity it is necessary for the weighting functions to be functions 
of {R t / R(x)} ^ 1 , with q > 1 [4], The weight functions are of the form 


W = Q(x)e imd 


R, 


R(x) 


g-i 


fiVW*) * <1*1 ) 


D(x)P(x)e ,mQ 


(17) 


where 


D(x) = 


R, V 

,R(x)j 


(18) 


In the parent element 


D{Z,,r\) 



(1 -t) q 


(19) 
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The weight functions are the complex conjugates of the shape functions multiplied by the 
additional decay term. 

THE WEAK FORMULATION IN THE INFINITE ELEMENT REGION 

The weak formulation of equation (2) for the infinite element region in which the steady 
flow is necessarily uniform is obtained by using equations (12) and (17) defining the assumed form 
of solution and the weight functions in the infinite element region. The gradient operations on the 
assumed shape and weighting functions yield 

V<|) = (V? - n\ r PV\i)e~‘^ (20) 


and 


1W = {D VP * + /r) r £)PVp. + PVD) e 




( 21 ) 


The superscript (*) denotes the complex conjugate of the operation 
\jp - p { + p e - i — P e Q , which is defined explicitly for the cylindrical coordinate system and 
takes into account the angular harmonics depending on m By using standard finite element 
operations equation (2) can be formulated to yield complex element “stiffness’ matrices K j} 
defined in terms of real mass, stiffness and damping matrices 


K 




( 22 ) 


where 


dP dP dP. dP, m 2 

V 




+ p ,[( 1 


dx dx dr ar 


(23) 
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(24) 


U, - [ff D r,r, ( 1 - [( 1 - M>) (|r- * <f )* * 2A/| ] } dtr 


c ,-ff now-** 


fl[(1 

[(1 - M 2 ) 


dx 3x 3x dr dr * J 
dD + M— + — P P } dV 


Sx ox 


0x 3r 3r 


(25) 


The definitions of the stiffness, mass, and damping matrices of equations (23)-(25) are 
implemented at the element level using the infinite mapping to the parent element. These results 
reduce to those of reference [4] when the medium is stationary and when account is taken of the 
operations which are particular to the cylindrical coordinate system. 

The mass matrix of equation (24) vanishes if the surface C r separating the standard finite 
element region from the infinite element region is a surface of constant phase for an apparent 
acoustic source location x 0 , r Q which is common for all elements. This is shown by referring to 
the definition of p(£), 


p(x) = Tl r (l|f(x) - ^l) 


(26) 


where 


i|r(x) = -1[-M(x -x 0 ) +R] 

P 2 


(27) 


and 


R = y/(x - X 0 ) 2 + p 2 (r - r 0 ) 2 (28) 

Since it is stipulated that C r is a constant phase surface, it follows that \jr 1 is constant. The 
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apparent source location is the same for all elements, leading to the conclusion that x 0 , r Q are 
constants. It can then be verified that 


(1 - m 2 )(^) 2 + (^) 2 +2 = 1 

dx dr ox 


(29) 


which from equation (24) leads to the result = 0 . This is consistent with the findings of 
Astley, Macaulay, Coyette, and Cremers [1] in the case of a stationary medium when the surface 
C r is a sphere, a constant phase surface in this case. While of some interest in the time harmonic 
formulation considered here, the vanishing of the mass matrix is of central importance when a 
time dependent formulation is implemented in the stationary medium case [2], It remains to be 
established that this is equally important in the case of a uniformly moving medium. It should be 
noted that the result M tj = 0 assumes that FEM interpolations and integrations are exact. There 
is in fact approximation error which is small as verified in calculations reported here. However in 
transient calculations this should be considered [5], 

TURBOFAN NACELLE EXAMPLES 

Computational examples will be given to demonstrate the performance of mapped infinite 
wave envelope elements as compared to standard wave envelope elements and in particular the 
performance of the new elements as a function of the expansion order will be examined. The 
codes which uses standard wave envelope elements [6,7] have a substantial history of 
benchmarking against experiment and simple test cases and can comfortably be used as a basis for 
evaluating the mapped elements 

Figure 1 shows the generic geometry of a turbofan inlet in an x , r plane of a cylindrical 
coordinate system. The acoustic source is on the plane Cy and produces a combination of radial 
modes at a specified angular mode m and non-dimensional frequency r) r . The source strength 
is specified by complex modal amplitudes [6], The nacelle has a forward velocity specified by the 
Mach number M 0 , which is represented for the stationary nacelle by a steady flow directed 
toward the nacelle. The steady flow into the nacelle is specified by the Mach number M i , taken 
to be uniform on the source plane. The steady flow field inside and outside the nacelle computed 
on the FEM acoustic mesh provides data for the FEM acoustic calculations. Figure 1 indicates 
that the computational domain is limited by an artificial baffle C b . Acoustic radiation is highly 
directional and it has been shown that the baffle can be oriented to have only minimal influence 
on the radiated field. This baffle is introduced to limit the dimensionality of the FEM 
discretization. 


10 


59 



The details of the FEM acoustic computations with the domain closed by a conventional 
wave envelope transition region to a Sommerfeld radiation boundary are given in [6], In a first 
example given here the propagation and radiation problem is formulated with a standard FEM 
discretization in the steady non-uniform flow near field and the domain is completed in the far field 
by the use of standard wave envelope elements. The specific case shown is at a reduced frequency 
T] = 25 and angular mode m = 23 with only the first radial mode incident. Only one radial mode 
propagates and it has a cutoff ratio near unity, which suggests that the peak lobe of the radiation 
pattern will be at a high angle relative to the nacelle axis. Steady flow is defined by M 0 = 0.3 and 
M = 0.2 . This case could correspond to rotor noise for a mildly supersonic tip speed rotor with 
23 blades. The conventional FEM mesh for this example is shown in Figure 3 (as used in the 
infinite element implementation). Infinite elements provide the reflection free boundary condition 
on the outer boundary. The domain limiting baffle is swept back more than 130 degrees from the 
duct axis. The mesh refinement is at close to the limit for the non-dimensional frequency 
considered. The proximity of the boundary C r to the nacelle is limited by the extent of the non- 
uniform flow field generated by flow around the nacelle. 

The results which are displayed are contours of constant acoustic potential magnitude in 
an x ,r plane superposed on the nacelle geometry. Figure 4 , for wave envelope closure of the 
computational domain, shows the radiation pattern in terms of contours of constant acoustic 
potential magnitude in the acoustic near field, that is, in the standard FEM region. Acoustic 
contours show a single lobe of radiation with the highest level on the transition surface C r 
corresponding to the closed contour. Considerable evidence of reflection from the wave envelope 
element region is revealed by the wavy quality of the iso-potential contours. In this case the 
transition surface C r (a constant phase surface) intersects the x axis at 2.5 duct radii.. This 
distance is apparently insufficient to achieve a good non-reflecting boundary using the wave 
envelope elements, which in this case consists of seven layers of elements extending to 10 duct 
radii on the x axis where the Sommerfeld condition is imposed. Mesh refinement has virtually 
no effect on the quality of the solution. The case depicted with radiation to the sideline is difficult 

because of the geometrical complexity of the source. 

Figure 5,6, and 7 show similar results when mapped infinite wave envelope elements of 
order 8 , 9, and 10 (the asymptotic expansion in powers of RJ R up to 8, 9, 10) are used to 
provide a reflection free boundary. In this case the transition boundary C T intersects the x axis 
at 2.0 duct radii, closer to the nacelle than for the results of Figure 4. The mesh refinement is 
approximately the same as used in Figure 4, and the element count is lower. There is 
progressively less evidence of reflection (waviness of the contours) as the element order is 
increased in successive figures. There is only a modest improvement between the order 9 
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expansion of Figure 6 and the order 10 expansion of Figure 7. The baffle limiting the 
computational domain, as noted in Figure 1, accounts for most of the residual evidence of 
reflection. Note that the iso-acoustic pressure contours of Figure 4 do not distribute in exactly the 
same way as those in Figures 5, 6, and 7 (they normalize on a different boundary). The result is 
that there are contours closer to the baffle in the latter cases, and these seem to show a small 
effect of the baffle. Timing in all of the cases shown is very nearly equivalent. 

The infinite element mapping and shape functions have also been used to generate the 
solution in the far field, and this is shown in Figure 8 where acoustic pressure contours are 
displayed. Acoustic pressure is obtained by post-processing acoustic potential [6.7] and is less 
accurate for a given mesh resolution because post-processing requires the spatial derivative of 
potential. In general, it is found that the solution for either potential or pressure is better in the 
infinite element region than in the conventional FEM region. This is because the infinite elements 
have shape functions which include the spatially harmonic character of the solution. The effect of 
reflection from the baffle is clearly shown, however the principal lobe of radiation is essentially 
unaffected by the baffle. 

As a final example to demonstrate the robustness of the reflection free boundary condition 
acoustic radiation from a turbofan exhaust is considered. Figure 9 shows the geometry with a 
potential flow jet issuing from the nacelle to represent exhaust flow. The mesh used in this 
implementation is shown in Figure 10. Note that there is a region of triangular elements that is 
used to work around the sharp tr ailin g edge of the duct without unacceptable element distortion. 
This has been found to be the most convenient mesh strategy consistent with the goal of keeping 
a simple, but optimal node numbering scheme for the frontal solver which is used. 

The FEM formulation in the interior region is somewhat more complicated than in the 
inlet case, requiring continuity of acoustic pressure and particle displacement across the shear 
layer separating the jet from the exterior flow [7], The length of the jet region forces the inner 
region in which standard FEM methods are used to be of much larger extent than in the inlet case. 
The case considered here is at non-dimensional frequency r\ r = 25 , with angular mode 
m = 23 , n = 1 incident. The jet Mach number is AT = 0.5 and the exterior flow Mach number 
is M a = 0.2 . The standard FEM region ends at 3.75 duct radii from the origin (on the axis of 
symmetry). The jet shear layer ends at 1 .0 duct radii from the axis origin, which corresponds with 
the tip of the center body. However, potential flow mixing persists for some considerable distance 
beyond this point. Figure 1 1 shows contours of constant acoustic potential in the near field 
generated with standard wave envelope elements providing the reflection free boundary 
condition. Figure 12 shows similar contours generated with the domain closed using tenth order 
mapped infinite wave envelope elements. The improvement with the mapped infinite elements is 
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substantial 

Some interesting properties of the solution are revealed in Figures 11 and 12. In the 
interior of the duct evidence is seen of standing waves along the duct wall. This is due to 
reflection at the duct termination. Discontinuity of acoustic potential (required by the physics of 
the problem), though not great, is seen across the shear layer between the jet and the surrounding 
flow. This is difficult to see due to the size of the figures, but it can be observed most clearly near 
the duct lip. The contours shown cover a range of 40 dB. There is a substantial improvement in 
the quality of the contours generated with the mapped elements. This is in spite of the fact that 
the mesh density used is marginal for the non-dimensional frequency considered and the fact that 
with the extended jet the boundary of the conventional FEM region is very close to the “extended” 
body. Figure 13 shows the post-processed acoustic pressure solution in both the near and far field. 
The effect of the baffle is noted to not substantially alter the principal lobe of radiation. 

CONCLUSION 

With suitable modifications mapped infinite wave envelope elements can be used to 
provide an effective reflection free boundary for acoustic radiation in a uniform steady flow. The 
fundamental solution for a source in uniform flow forms the basis for an asymptotic expansion in 
R-i in the infinite elements, where R is the “converted radius”, R 2 = x 2 + P 2 r 2 . The order 
of the asymptotic expansion can be chosen to meet the needs of the problem. Element mapping 
functions are identical to those previously proposed for the stationary medium case and the shape 
functions are of the same form as those in the stationary medium case with differences only in the 
details. Examples show that mapped infinite wave envelope elements provide a superior 
reflection free boundary for cases in which standard wave envelope elements generate reflections 
which appear in the radiated field. It has been demonstrated in the nacelle inlet case that this 
improved reflection free performance can be achieved on a reduced mesh. 
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C b is an artificial baffle to limit the computational domain with minimal 
influence on the solution. 
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Figure 3. Standard FEM mesh used for inlet radiation showing domain limiting baffle. Radiation 
condition on the outer boundary implemented with mapped infinite wave envelope 
elements. 
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Figure 9. Geometry of a generic exhaust duct showing the shear layer separating the exhaust 
flow and external flow. 
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Figure 10. Standard FEM mesh used for aft fan duct radiation showing domain limiting baffle 
and region of triangular elements. Mapped infinite wave envelope element closure on 
the outer boundary 
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Abstract 

An explicit, analytical, multiple-scales so- 
lution for modal sound transmission through 
slowly varying ducts with mean flow and acous- 
tic lining, is tested against a numerically “ex- 
act” finite element solution. The test geometry 
taken is representative of a high-bypass turbo 
fan aircraft engine, with typical Mach numbers 
of 0 . 5 - 0 . 7 , circumferential mode numbers m of 
10-40, dimensionless wave numbers of 10-50, 
and both hard and acoustically treated inlet 
walls of impedance Z = 2 — i. Of special inter- 
est is the presence of the spinner, which incor- 
porates a geometrical complexity which could 
previously only be handled by fully numerical 
solutions. The results in predicted power atten- 
uation loss show in general a very good agree- 
ment. The results in iso-pressure contour plots 
show good comparison in the cases where scat- 
tering into many higher radial modes can occur 
easily (high frequency, low angular mode), and 
again a very good agreement in the other cases. 


Introduction 

The calculational complexities of the multi- 
ple-scales solution for modal sound transmis- 
sion through slowly varying ducts with mean 
flow and acoustic lining (presented in [ 1 ]), are 
no more than for the classical modal solution 
for a straight duct. The multiple-scales so- 
lution is an approximation utilizing the axial 
slope of the duct walls as small parameter. This 
slope is for aerodynamical reasons indeed in- 
variably small in any aero-engine duct. 

Therefore, this multiple-scales solution pro- 
vides an interesting alternative in aero-engine 
applications, as it both allows for the advan- 


tages of the analytical approach (speed of calcu- 
lation and relative simplicity of programming), 
and variable geometries including spinner and 
mean flow variation. 

The final approximation error in realistic 
geometries, however, is difficult to determine, 
except for an order of magnitude estimate say- 
ing that it scales on this slope. It is therefore 
of interest to directly compare the analytical 
approximation with a state of the art fully nu- 
merical solution of the same physical model. 
This is the subject of the present paper. 

As a first step towards exploring the possi- 
bilities, a series of tests are carried out, com- 
paring the analytical results with results of the 
finite element solution, given in [ 2 ], of a com- 
pressible inviscid isentropic irrotational mean 
flow, superimposed by linear acoustic pertur- 
bations. 

Physical model 

We consider a circular symmetrical duct 
with a compressible inviscid perfect isentropic 
irrotational gas flow, consisting of a mean flow 
and acoustic perturbations. To the mean flow 
the duct is hard-walled, but for the acoustic 
field the duct is lined with an impedance wall. 
In view of the adopted aero-engine geometry, 
the inner wall (corresponding to the spinner) 
will be hardwalled, without lining. 

We make dimensionless: spatial dimensions 
on a typical duct radius Rqq , densities on a ref- 
erence value pco t velocities on a reference sound 
speed Coo, time on Roo/cco , pressure on 
and velocity potential on RoqCcg • Note that 
the corresponding reference pressure poo satis- 
fies pcoC^ = 7 Poo, where 7 = 1.4 is the (con- 
stant) ratio of specific heats at constant pres- 
sure and volume. 
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The fluid in the duct is described by 

Pt + V* (pv) =0, (la) 

p(v t + v- Vv) + Vp = 0, (lb) 

7 p = p 7 , c 2 = = P 7-1 (lc,d) 

(with boundary and initial conditions), where 
v is particle velocity, p is density, p is pressure, 
c is sound speed (all dimensionless). 

Since we assumed the flow to be irrotation- 
al, we may introduce a velocity potential 0, 
such that v = V0, and the above momen- 
tum equation may be integrated to a variant 
of Bernoulli’s equation 


impedance Z. The pertaining boundary condi- 
tion is for a point near the wall but still (just) 
inside the mean flow. For arbitrary mean flow 
along a (smoothly) curved wall, with normal n 
directed into the wall, this was given by Myers 
[3], eq. 15, as 

iu;(v-n) = j^iw 4* V* V — n-(n* VV)j • (5) 
Geometry 

The reference values taken for non-dimen- 
sionalization are at the source plane x = 0, 
including the outer radius for length scale. 

The outer radius Rn and inner radius R\ 
are described by the following formulas 


<90 

~dt 


+ 2 |V| ' + 



= const. 


( 2 ) 


This flow is split up into a stationary mean 
flow part, and an acoustic perturbation. This 
acoustic part varies harmonically in time with 
circular frequency u;, and with small amplitude 
to allow linearization. 

In the usual complex notation we write then 


R 2 (x) = 1 - 0.18453x /2 

e -n(i-V) _e“ n 

4-0.10158 n , (6) 

1 — e _li 

Ri(x) = max[0, 0.64212 

-(0.04777 + 0.98234r ,2 ) 1/2 ], (7) 

where x' = xjL and L = 1.86393 is the length 
of the duct; see figure (1) 


v = V + ve'"*, 4> = $ + <£e ,u,t , 
p=D + pe iut , p = P + pe lult , 
c= C + ce 1 ^. 

Substitution and linearization yields: 


• mean flow field 

V.(DV) = 0, (3a) 

ll v | 2 + — T = E ( 3b ) 

2 7—1 

C 2 = 1 P/D = D 7_1 ; (3c) 

• acoustic field 

iujp 4- V • (DV0 -f pV) = 0, (4a) 

iw<?> + V-V^.+ ^ =0, (4b) 

p = C 2 p, (4c) 



where E is a constant, and the acoustic pertur- 
bation of c is further ignored. The integration 
constant in equation (4b) may be absorbed by 
0. For the mean flow the duct wall is solid, so 
the normal velocity vanishes at the wall. The 
subsonic mean flow is determined by conditions 
of uniformity upstream, the constant E , and an 
axial mass flux 7 tF . 

For the acoustic part the outer duct wall is 
a locally reacting impedance wall with complex 


Figure 1: Geometry 

The mean flow is selected such that at the 
source plane x = 0 the Mach number is equal 
to —0.5, and the dimensionless density equal to 
1. The corresponding axial Mach number and 
dimensionless density variation (based on the 
quasi-one dimensional mean flow solution; see 
below) is depicted in figure (2). 
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Figure 2: Mach number and density 


Multiple scales solution 

For the success of the analytical solution it 
is essential that the mean flow and the acoustic 
field are approximated on the same footing. An 
arbitrary, ad-hoc, mean flow field would not 
allow the present explicit solution. So the mean 
flow used for the multiple scales solution is not 
exactly the same as the one used for the finite 
element solution. They are, however, in terms 
of approximation of the same level. Therefore, 
we give here the mean flow and the acoustic 
field together. 

The approximation is based on the assump- 
tion that geometry and mean flow vary slowly, 
he. on a length scale much larger than a typ- 
ical duct diameter or wave length. This is, of 
course, for aerodynamical reasons the case in- 
side an aero-engine inlet duct. We introduce 
the ratio between a typical diameter and this 
length scale as the small parameter 5, and re- 
write the duct surface (in radial coordinates 
(x, r, #)) 

r — Ri(X), r=R 2 (X) } X = ex. (8) 

By rewriting Rx i2 as a function of slow variable 
X , rather than x, we have made our formal 
assumption of slow variation explicit in a con- 
venient and simple way. Although in the final 
result 6 will play no explicit role, a represen- 
tative value of e will be necessary for an order 
of magnitude estimate of the approximation er- 
ror. 

By assuming that the mean flow is nearly 
uniform with axial variations in X only, we find 
that small axial mass variations can only be 
balanced by a small radial flow, so 

V - U 0 {X)e x +eVi(X i r)e r 


and similarly are P ~ Pq(X), D ~ Dq( A), and 
C ~ C 0 (X) to leading order only dependent on 
A'. It follows that 

F 

Uo{X) ~ Do(X)(R 2 2 (X) - R\(X)) 

with V i , D 0 , Pq and C 0 are given by the well 
known one dimensional gas flow equations (see 
e.g. [1]). t 

The acoustic field is assumed to be described 
by mode-like solutions of the form 

4>(x, r, 9; e) = .4(A\ r; e) 

( 9 ) 

After expanding A = Aq + eAi + 0{e “) and /j = 
^ Q _l -0{e 2 ) (any possible pi can be absorbed by 
Ai), and substitution in equations and bound- 
ary conditions, we find for Ao a Bessel-type 
equation in r, so we obtain the slowly varying 
mode 

Aq = N{X)J m {Q{X)r) + M(X)Y m (a(X)r) 

( 10 ) 

where J m and Y m are the m-th order Bessel 
function of the first and second kind. Radial 
eigenvalue a and M/N are determined by the 
boundary condition, while 

a 2 + jj? — Dr /Cl , Q = id — fj.Uo. 

The crux of the solution is the determination 
of amplitude N(X), as a function of X . This 
is determined by the next order equation for 
A\. It is, however, not necessary to solve this 
complicated equation. A solvability condition 
[1] is enough to generate a differential equation 
in X for N, which can be solved exactly. The 
general solution for the hollow cylinder (R\ = 
0, M = 0) is given by 


(Qo \ 2 _ 

( DoujaRl f 

O -0 

m- - \ 

V N ) 

{ 2 Co i 1 

a-R\ > 


+ ^ C 2 V m(o/i2)2 (n) 

where Qq is an integration constant, and 
f 2 = Q 2 DoRi/iuZo, cr 2 = l-(C 0 2 -^ 0 2 )a 2 /u; 2 . 
The solution for the annular cylinder is more 
complicated, although explicit, and can be 
found in [1] 

Finite element solution 

A numerical model for sound propagation 
is based on a finite element discretization of 
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the steady flow field equations (3a-3c) and the 
acoustic field equations (4a-4c). The weak for- 
mulation of equation (3a) for the steady com- 
pressible flow in the duct is in terms of the 
steady flow velocity potential and steady flow 
density D , 



VW-(DW) dV 



s 


W(DV<f>)-ndS . 

( 12 ) 



Figure 3: Computational domain FEM 

Weighting functions W are from the class 
of continuous functions on the volume V of the 
duct bounded by the duct surface 5, which 
includes the duct wails and source and exit 
planes. A solution for $ is sought in the class 
of continuous functions. The unit normal is 
directed out of the duct. The duct geometry 
and the steady flow field are axially symmetric 
favoring the introduction of a cylindrical coor- 
dinate system with x axis coincident with the 
axis of symmetry, r axis in the source plane 
at x — 0, and the angular coordinate 9 locat- 
ing the r axis in the x — 0 plane (see figure 
3 for the computational domain). The steady 
flow field is represented in an (z, r)-plane, and 
is two-dimensional. A standard finite element 
formulation of equation (12) is based on eight 
node isoparametric serendipity elements. 

Equations (3b) and (3c) are subsidiary re- 
lations that are used in an iterative solution 
in which at each stage the finite element dis- 
cretization of equation (12) is solved with a 
density and speed of sound field derived from 
the previous iteration step. The boundary in- 
tegral on the right hand side, which is a natu- 
ral boundary condition, specifies the mass flow 
rate on the source plane. A forced boundary 
condition setting the level of the potential is 
required on the exit plane (figure 3). It is as- 
sumed that the source plane and exit plane are 
located remotely enough from regions of non- 
uniformity in the duct so that at the source and 
exit planes the flow velocity is uniform, per- 


mitting the natural and forced boundary con- 
ditions to be easily implemented. In practical 
calculations in the type of duct considered this 
turns out not to be restrictive. 

The mean flow' in the duct given by figure 
3, with a uniform Mach number at the source 
plane M = -*0.5, is directed from right to left 
(an inlet flow). The duct shape defined by 
equations (6) and (7) begins at the source plane 
x = 0 and is extended beyond the nominal ter- 
mination used in the analytical development in 
a uniform duct to allow the flow field to become 
uniform at the exit plane. No extension is used 
at the source end and the extension at the exit 
end is probably longer than necessary. V$, C, 
D are required data for the FEM solution for 
acoustic propagation. 

A finite element model for acoustic propa- 
gation is based on a weak formulation of equa- 
tions (4a-4c). Acoustic perturbations in pres- 
sure, density and velocity potential are har- 
monic in time with frequency u and harmonic 
in the angular coordinate 9 of the form p(x } r), 
p(x, r), and <p(x, r) times the complex exponent 
e iwf-im0 The weak formulation [2] is 


III t vvy '( DV ^ + _ iuW p) dv 

= JJ W{DVJ> + pV<t>)-ndS (13) 


The weighting functions are taken as 
W[x,r)e' md . Angular harmonics proportion- 
al to e~ [me represent the decomposition of the 
solution periodic in 9 in a Fourier Series. The 
angular mode number m is a parameter of the 
solution. The surface integral is over all sur- 
faces bounding the domain. The unit normal 
for the surface integral is out of the domain at 
the surface in question. The w'eak formulation 
continues with the linearized momentum equa- 
tion (4b) and linearized equation of state (4c), 
to rewrite equation (13) in the form 


III |H c2w ' v ^ _(v ' vw0(v ' v, ^ + 

V iw[w(V.V<*)-(V-VW>] -u> 2 W<j>} dV 

= II -^{c 2 WV<t>-VW(V-V<j>) 

—[uiVWip j -ndS (14) 


Note that the local steady flow dimensionless 
velocity V is equivalent to the reference Mach 
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number M r which in fact is local steady flow 
velocity divided by the speed of sound at the 
source plane. 

The surface integral on the right hand side 
of equation (14) is the natural boundary con- 
dition. On the duct walls this provides the 
boundary condition for either rigid walls (the 
integral vanishes) or for a normally reacting lin- 
ing with an impedance specified by equation 
(5). In the present FEM implementation equa- 
tion (5) is simplified by the elimination of the 
term involving n*(n*VV) on the right hand 
side. In the duct geometry studied here this 
term, which depends on the nonuniformity of 
the steady flow field, is of little importance in 
affecting attenuation (even though it is asymp- 
totically small but not negligible, and crucial 
in completing the analytical formulation). 

Details of the FEM procedure for discretiza- 
tion of equation (14), and other references, can 
be found in [2]. A discussion of the implemen- 
tation of the impedance boundary condition is 
available in reference [5]. 

On the source plane and exit plane the nat- 
ural boundary condition is used to introduce 
the noise source and non-reflecting boundary 
conditions. On these planes the acoustic poten- 
tial is recast via an eigenmode expansion such 
that the acoustic potential is given in terms of 
the complex amplitudes of the right and left 
propagating acoustic duct modes appropriate 
for the geometry and flow conditions which pre- 
vail there [6]. On the source plane, x = 0 in the 
present study, right propagating modal ampli- 
tudes at the source plane are specified via a 
forced boundary condition. Left running (re- 
flected) modal amplitudes at the source plane 
and right running modal amplitudes at the ex- 
it plane are unknown and part of the solu- 
tion. Left running modal amplitudes at the 
exit plane are forced to vanish, imposing a non- 
reflecting boundary condition. Details of the 
modal boundary condition are available in ref- 
erence [4]. 

Finite element discretization for acoustic 
propagation is carried out on the same grid 
with the same element type as used in the steady 
flow model. Required data generated in the 
steady flow representation is transferred direct- 
ly to the acoustic analysis. Mesh density is gov- 
erned by the demands of the acoustic problem 
and is substantially more refined than would be 
required in the steady flow analysis. 

The FEM solution proceeds with the com- 


putation of the acoustic potential field. Post- 
processing by the use of equation (4b) gener- 
ates the acoustic pressure field. The solution 
also includes reflected modal amplitudes and 
transmitted modal amplitudes. Acoustic power 
reflection and transmission characteristics are 
computed directly from the input modal ampli- 
tudes and computed reflected and transmitted 
modal amplitudes. Reciprocity characteristics 
of the scattering matrices and acoustic power 
balances are also monitored as a check of com- 
putational accuracy in the case of no acoustic 
treatment on the duct walls [7]. 

Post processed acoustic pressures are repre- 
sented on iso-pressure amplitude contour plots 
superimposed on the duct geometry. Compari- 
son of FEM and multiple scales results is based 
on visual comparison of these contours, but 
perhaps more importantly on the basis of com- 
puted power transmission coefficients. 

Differences between Multiple Scales 
and FEM Formulations 

There are minor differences between the 
multiple scales solution and the finite element 
model. The field equations (3a)-(3c) and (4a)- 
(4c) are exactly the same in both cases, includ- 
ing the convention for non-dimensionalization. 
The implementation of the impedance bound- 
ary condition in the FEM formulation neglects 
the n* (n- VV)-term in equation (5). This is 
done principally to simplify the FEM imple- 
mentation of this difficult boundary term (it 
requires a gradient of the steady flow velocity 
which is already the gradient of the steady flow 
potential). For a cylindrical duct this term is to 
leading order in e equal to ±V ~ eURx/R- So 
it is small, but asymptotically not smaller than 
any other effect due to duct variation. Never- 
theless, we found as yet no indication that the 
effect on attenuation predictions is significant. 

The FEM formulation includes the propa- 
gation of many modes and therefore scattering 
is an integral part of the solution. This is man- 
ifested by reflection of the incident mode and 
other modes which are not incident as well as 
the transmission of modes which are not inci- 
dent. The multiple scales method utilizes the 
fact that in the smooth parts of the duct any 
scattering into other modes is normally neg- 
ligible. The propagating sound is still mode- 
like, albeit not in the strict sense of self sim- 
ilar straight duct modes, but mode-like solu- 
tions with slowly varying amplitude and phase. 


5 


American Institute of Aeronautics and Astronautics 


81 


At abrupt changes in geometry scattering into 
other radial modes may be included (not done 
here) by methods like mode-matching. 

The FEM solution requires the source to be 
represented in terms of input modal amplitudes 
for eigenmodes for a duct with hard walls. The 
source is always located in a section of the duct 
which has rigid walls. The net effect is that 
there is always a transition from a rigid wall 
to an impedance wall at both ends of the duct 
in the FEM model. This has implications for 
scattering which are not readily quantifiable. 

In the multiple-scales analysis the gener- 
al solution is built up from a summation over 
slowly varying modes. The natural way to test 
its validity is therefore to study a single, soft- 
wall, mode. In order to generate an equivalent 
input in the FEM model it is necessary to rep- 
resent the soft wall eigenmode as an eigenmode 
expansion of hard wall eigen modes. Since this 
sound field distribution, presented to the lined 
duct, essentially “fits” directly into one soft- 
wall mode, only little reflection at the source 
plane is to be expected. A single mode multiple- 
scales analysis is therefore simulated by a mul- 
tiple mode FEM solution. 

Finally, it is noted that the FEM model re- 
quires conditions at the source and exit planes 
which give rise to reflected modal amplitudes 
and to the reflection free termination (or speci- 
fied or computed reflection characteristics). 
These conditions do not play an immediate role 
in the analytical solution, with left and right 
running waves already given explicitly. Al- 
though inherent in any practical application, 
we have not tried to model these conditions in 
the analytical part of the present tests in or- 
der not to obscure the comparison and to re- 
strict the sound field to that of a single, right- 
running, mode. 

It is not possible to make FEM and multiple- 
scales models exactly equivalent, nor should 
it be, since the multiple-scales solution is an 
approximation based on well documented as- 
sumptions. It is a goal of the numerical com- 
parisons to be given here to investigate how 
successfully the multiple-scales solution repre- 
sents the more exact FEM model. 

Results 

The cases considered are grouped as the fol- 
lowing 4 series of iso-pressure contour plots: 
the first radial mode of 


fig. 4 m = 10, u = 10, 

fig. 5 m — 10, oj = 16, 

fig . 6 m =: 10, uj = 50, 

fig. 7 m — 40, cj = 50, 

is input under the following conditions: 

fig. a hard wall, no flow, 

fig. b hard wall, flow, 

fig. c soft wall, no flow, 

fig. d soft wall, flow. 

The left column of the figures is the numerical 
FEM solution, the right column the analyti- 
cal MS (multiple scales) solution. “Soft wall” 
denotes a wall impedance Z = 2 — i. “Flow” 
denotes a mean flow with Mach number -0.5 at 
the source plane. “First radial mode” denotes 
in general the mode with smallest real part of 
radial eigenvalue q. For the soft walls the pre- 
dicted attenuation (10 log of ratio of acoustic 
power through source and inlet plane) is giv- 
en in the caption of the figures. For the hard 
walls the attenuation is either zero (mode is cut 
on) or infinite (mode is cut-off and reflection is 
negligible). 

The selected cases do not show turning point 
behaviour (hard-wall cut-on, cut-off transition). 

The possible differences between FEM and 
MS are due to the following errors or modelling 
discrepancies. 

1. The approximation error of 0(e 2 ). For this 

we need an estimate of e. Suitable is a typical 
value of ~R 2 = = 0(e). From for- 

mula (6) it appears that R!n varies between 
-0.12 and 0.12 along [0,1.75], but increases 
to 0.4 in the lip region [1.75,1]. If we take 
e =: 0.1, the estimated approximation error 
is a few percent. 

2. Small but inherent reflection in FEM at inlet 
plane and lip region. 

3. Not exactly the same source in soft wall cas- 
es, since FEM uses a source defined by an ex- 
pansion in a finite number (15) of hard wall 
modes, a small distance ( T3o avva y ^ rom 
the lined section. 

4. Slightly different impedance definition in flow 
cases. 

5. Slightly different mean flow. In MS the mean 
flow field is approximated to the same level 
as the acoustic field. 

For highly attenuated or only cut-off modes 
(m = 10, uj = 10) of fig. 4a-d. the agreement is 
almost perfect. None of the above errors seem 
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to play a role. 

For the series (m = 10, u; = 16) of fig. 5, 
with 1 (no flow) or 2 (flow) modes cut-on, the 
agreement is good in the iso-contour plots, and 
almost perfect in attenuation. Some wiggles 
are visible in fig. 5b what is probably error type 
1 T due to interference with the other cut-on ra- 
dial modes. 

The high frequency series (m = 10, u = 50) 
of fig. 6 has very low acoustic pressure val- 
ues near the duct axis, and many radial modes 
cut on (at source plane, no flow: 9, flow: 11). 
We see strong interference with these higher 
modes. The most important region near the 
outer wall, however, is in very good agreement, 
and for no flow the attenuation agrees exact- 
ly. With flow the 1st and 2nd radial soft wall 
modes happen to be rather close to each oth- 
er, and a residual second mode (due to error 
3) leads to a slightly (0.5 dB) different attenu- 
ation. 

The high m, high frequency [m — 40, u; = 
50) series of fig. 7 has 2 (no flow) and 3 (flow 7 ) 
radial modes cut on, giving rise to some wiggles 
in figs. 7a, b. Effects due to error 3 are probably 
visible in fig. 7d, although the predicted attenu- 
ation agrees very well. The (academic!) differ- 
ence in attenuation of fig. 7c (195 and 210 dB) 
is no numerical round-off error, but due to the 
fact that the plotted mode is least attenuated 
at the inlet but not at the source plane. Residu- 
al modes due to error 3 are likely to dominate in 
the FEM solution near the source plane, lead- 
ing to a different attenuation. 

Conclusions 

Any selection of test cases is necessarily lim- 
ited. It would have been easy to create a more 
or a less favourable comparison, by making a 
suitable selection of geometry and parameters. 
This is not done here. We have defined the test 
runs entirely on the basis of their relevance to 
turbo fan engine inlet duct applications, and 
we have not skipped unfavourable cases after- 
wards. The only restriction we made was that 
no cut-on/off transition in the hard walled duct 
be present. This phenomenon is not yet includ- 
ed in the analytical solution, while at the same 
time, of course, it is absent in any lined duct. 

So considering the fact that the cases are 
likely to be a representative cross section of re- 
ality, we think the conclusion is justified that 
the MS and FEM solutions compare favourably. 


both in iso-pressure contours and in predicted 
attenuation. Principle differences are related 
to scattering at inlet plane, and to input mode 
synthesis. The best results are obtained with 
lining (reducing importance of reflection) and 
when the modal structure permits few or no 
cut-on scattered modes. The attenuation dif- 
fers in general no more than a few tenths of a 
dB. 

The correlation shows that MS is definitely 
useful in applications for assessing liner perfor- 
mance in realistic geometries. Both extending 
the theory, and further comparison with FEM, 
for example with an MS implementation that 
includes a complete modal spectrum and open 
end reflection, is therefore to be recommended. 
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ABSTRACT 

In order to achieve satisfactory results 
with finite element formulations for pressure in a 
propagating acoustic field with relatively high 
non-dimensional frequencies, the problem of 
ever-diminishing mesh size must be resolved. 
The convected potential formulation describing 
an acoustic field introduces problems not 
addressed in the simple Helmholtz equation. 
Post-processing to calculate pressure is necessary 
and this leads to additional dimensionality 
problems beyond those encountered in modeling 
acoustic potential due to an inaccuracy appearing 
in the calculation of the potential derivatives. 
This oscillating erroneous behavior is rooted in 
the element shape functions and modifications 
have been made, using elements of higher order, 
to contain this discrepancy enough to where 
post-processing does not add significantly to the 
dimensionality problem. In so doing, satisfactory 
pressure models for non-dimensional frequencies 
up to 100 in a variable area circular duct with a 
wide range of subsonic Mach numbers and 
angular modes can be calculated in a reasonable 
sized domain. 

INTRODUCTION 

Application of the finite element 
method (FEM) in modeling acoustic propagation 
leads to problems of large dimensionalities when 
high frequencies are considered. A mesh of 
about 10 nodes per wavelength is considered 
adequate resolution in order to achieve 
acceptable results in modeling propagating 
waves. Simply reducing the mesh size is 
certainly an option. However, it becomes quickly 
impractical and computationally expensive. 
Several approaches have been studied in order to 
deal with the dimensionality problem. 

One common approach is to lower the 
number of nodes required per wavelength. In 
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acoustic radiation to the far field the 
development of infinite elements, which 
incorporates spatially harmonic radiation mto the 
shape functions, has greatly reduced the far field 
dimensionality problem' 0 . The earliest approach 
was based on exponential decays with radial 
distance 1 . Later, and with more success, mapped 
infinite elements were developed based on the 
asymptotic behavior of the far field and the 
correct spherical or cylindrical decay was 
incorporated 2 . Bumetri extended this method to 
use prolate spheroidal coordinates, allowing it to 
be more flexible in modeling objects with large 
aspect ratios, i.e. reducing the extent of the 
conventional mesh in the near field, thus 
reducing the overall mesh. A variation of the 
infinite element method is that of the wave 
envelope element method which restricts the 
computations to a large but finite domain . 
Astley 4 incorporated the use of complex 
conjugates weighting as opposed to Galerkin 
weighting, reducing the integrals to simpler 
forms, thus enabling Gauss-Legendre 
integration. The penalty was a non- symmetric 
coefficient matrix, although with less frequency 
dependence. This method was also capable of 
reducing the reach of the conventional element 
domain while extending into flow problems . 
Mapped infinite wave envelope elements have 
attributes of both infinite and wave envelope 
elements. This technique has recently been 
extended to include three-dimensional elements 
of variable order 3 , sperhoidal elements of 
variable order to improve modeling around a 
slender or flat object 9 , and a non-reflective 
boundary with uniform steady flow . The 
mapped infinite wave envelope elements can be 
used in what is normally considered the acoustic 
near field because of their adjustable 
interpolation order, thus reducing mesh 
refinement and dimensionality to a greater extent 
than with previous approaches. Chadwick and 
Bettes 11 modeled the phase p and the wave 
envelope A rather than the potential <|> (the 
potential of a traveling wave is expressed as 
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(|)=Ae‘ p ) for finite wave envelope elements and 
were able to use a coarser mesh in the near field. 
This was later extended to infinite elements \ 
The penalty however is that an iteration 
procedure is necessary. All of the previously 
mentioned approaches are well suited in an 
extenor region. Inside a duct, however, the 
complexity of the acoustic field makes for 
limited success. In fact, Astley and Eversman' J 
suggested the use of traditional wave envelope 
elements inside a duct before they were used in 
the far field, but it was effective in only special 
circumstances. 

The most traditional FEM approach in 
dealing with short waves is to let the mesh size 
decrease while the piecewise polynomials of 
degree p are constant. This is known as the h- 
version, while the p-version leaves a fixed mesh 
and allows p to increase. These have been 
combined and extended into the adaptive hp- 
version. For a detailed background and method 
of implementation see references [14] and [15]. 
In the field of acoustics this method has recently 
been used in modeling the Helmholtz equation 
for non-dimensional frequencies up to 30 tu 16 . 
Although a significant decrease in total degrees 
of freedom is achievable because of its adaptive 
behavior, it is computationally expensive due to 
the iteration procedure, despite exponential 
convergence. The nodal density is still needed at 
about 10 nodes per wavelength for reasonable 
errors. 

The incorporation of wavelike behavior 
into the elements has not been limited to infinite 
and wave-envelope elements. Babuska and 
Melenk introduced the partition of unity method 
(PUM) based on introducing prior knowledge 
about the governing differential equation into 
non-polynomial functions used for the solution 
approximation, rather than traditional (mapped) 
polynomials 1 '. These special functions 
approximate well the exact solution and the wave 
direction at each node is implicitly determined. It 
is advantageous to use a set of plane wave 
solutions of the homogenized Helmholtz 
equation as the local function basis. 
Computational cost is still a problem compared 
to non-iterative approaches. In reference [17] it 
is also suggested that PUM is suited well for a 
mesh-less formulation. More recently, an 
element-free one-dimensional Galerkin method 
was applied to the Helmholtz equation 18 . 

The combination of PUM and a 
standard FEM has been called the generalized 
finite element method and seems to offer a 
significant reduction in dimensionality versus 


standard FEM without requiring re-meshing 19 . In 
a recent presentation Bettes and Laghrouche 
extended this idea to the Helmholtz short wave 
problem 20 . An approach based on a somewhat 
similar idea as that behind the PUM was that of 
residual-free bubble functions -1 . It added these 
functions to the piecewise linear polynomials 
and used subspaces where these functions 
satisfied the differential equation strongly and 
solved analytically for them. This was applied to 
the Helmholtz equation for relatively low non- 
dimensional frequencies (<10). Franca and 
Macedo extended this to a more flexible two 
level method 22 . They used a submesh defined in 
each element interior and solved the differential 
equation numerically for the bubble functions, 
allowing for irregular meshes. 

It is important to note that the methods 
used in reference [14] to [22] have all been 
applied when solving the Helmholtz equation. 
The introduction of a moving medium in the 
acoustic field has lead to a different approach. 
The FEM model used in this paper is potential 
based and post-processes to find the acoustig 
pressure. This convected potential formulation 
introduces problems not encountered in working 
with the Helmoltz equation. Among them is the 
consideration of an effective wavelength, due to 
flow velocity and while the potential solution 
abide by the expected nodes per wavelength ratio 
of 10 to 12, the calculation of pressure which 
involves the derivative of the potential solution 
apparently does not. The nodal density required 
by the post-processing adds greatly to the 
dimensionality of the problem. An investigation 
into the cause of the shortcoming and possible 
simple modifications for this model was 
conducted. The variation in the order of the 
polynomials was performed for both a quasi-one 
dimensional and an axially symmetric system. 
The results were considered in a ducted acoustic 
field with non-dimensional frequencies up to 100 
and a wide range of propagating modes. It 
involved large area variations and most of the 
subsonic flow range. The behavior of this high 
frequency complex acoustic field was studied 
and the post-processing problem was greatly 
diminished while a possible method for an 
extension to improve the dimensionality of the 
overall system was considered. 

FORMULATQN OF PROBLEM 
The problem in question is the 
propagating acoustic field in a moving medium 
inside a circular non-uniform duct. A typical 
geometry of such a duct is seen in figure 1. The 
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noise source is at the left end and the medium 
may flow in either direction with left to right 
considered positive direction. The acoustic field 
description is based on a convective potential 
formulation obtained by considering unsteady 
acoustic perturbations on a steady compressible 
potential flow. Eversman et. al. previously 
developed a FEM model for the turbofan engine 
inlet radiation problem 5 ' 5 . This FEM model is 
basically the same without the exterior region but 
with the capability of handling compressible 
flow. An extensive description of the 
formulations and the computational scheme can 
be found in reference [6] and the references 
therein. However, a summary of the major 
features of the approach will follow. 



Figure 1. Non-dimensional computational 
domain 


It is assumed that the moving medium is 
non-viscous and that all processes are isentropic. 
The field equations for both the potential flow 
and the acoustic perturbations in non- 
dimensional forms are derived from the non- 
dimensional continuity, momentum, and energy 
equations, the latter in the form of the isentropic 
equation of state: 

^. + V(pK) = 0 (1) 

dt 

T_ + (K-V)F = V p (2) 

dt p 

P = V 0) 

7 

p , p , and V are non-dimensional pressure, 
density, and velocity, respectively, y is the ratio 
of specific heats. The non-dimensional speed of 
sound is 


P 

p r is the reference density and c r is the 
reference speed of sound and are taken as their 
respective values at the source plane. The 
reference length R is defined as the duct radius 
at this plane. Pressure has been made non- 

dimensional by p r c* , density by p r , velocities 
and speed of sound by c r , velocity potential by 
c r R , and time by Ric r . The introduction of 


velocity potential and the linearization of the 
conservation equations to the first order in 
acoustic perturbation yields 

V(pV^ 0 ) = 0 (5) 


for steady mean flow, and 

— + V-(A>V* + pV* > ) = 0 (6) 

dt 

for the acoustic perturbations, where <p Q and Po 
are the mean flow potential and density, 
respectively, and <j) and p are acoustic 
potential and density. The mean flow and 
acoustic densities are determined by 


Po = 


1 - V0 O ■ V fa) 


(7) 


P = ~ 


Po 

Co 2 


4r + (V*> •**) 

at 


( 8 ) 


where c$ = p$ 1 is the local speed of sound. The 
acoustic pressure is related to the acoustic 
potential by 




-Po 


+ ( V A) * ^ 0) 

dt 


(9) 


A standard finite element Galerkin 
approximation is used to solve the problem. The 
computational domain is modeled using either 
quasi-one-dimensional or axially symmetric 
formulations with mesh elements of different 
order. The solution is achieved through three 


steps: 

• The time invariant mean compressible flow 
problem. Equation (5) is solved in the 
domain in an iterative procedure for the flow 
potential. 

• Uniform duct eigenvalue problem. A 
uniform duct eigenvalue problem is solved 
on the fan face in order to express the fan 
and exit face boundary conditions in terms 
of duct eigenfunctions. 

• Acoustic propagation problem. Equation (6) 
is solved for the acoustic potential, using 
equation (8). The solution is desired in the 
case of a harmonic source on the source 
plane with time and angular dependence 

given by e i(n,r ' md) where rj r is the non- 
dimensional input frequency ( rj r =coR/c r , 
co is the input frequency) and m is the 
angular mode number. 

The acoustic pressure distribution in the duct is 
calculated through post processing of the 
acoustic potential solution, using equation (9). 
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Figure 2a. Pressure contours using quadratic elements tor q r — 25, M[ .l,m 1, and N/X 11.1 



Figure 2b. Pressure contours using cubic elements for rj r - 25, M\ - .1, m 1, and N/X 1 1.3 

moving medium. This is the source ot the 


As far as dimensionality is concerned the 
difficulties of this model lie withm this process. 
For high values of input frequency the nodal 
density of the mesh necessary to achieve 
acceptable pressure results has to significantly 
exceed 10 nodes per wavelength. An attempt to 
reduce the overall dimensionality problem can 
not be made until this ratio applies to satisfactory 
pressure calculations. 

RESULTS AND DISCUSSION 

In this section some numerical results 
are given to demonstrate the improvements made 
in reducing dimensionality of the domain. The 
solution of the acoustic radiation problem yields 
the acoustic velocity potential at the nodes of the 
mesh. The solution that yields the acoustic 
pressure is obtained from equation (9) as 

P = ~Pol ir lr& + ( V ^0 ' ( 10 ) 

Pressure at a node is calculated from the 
elements sharing that node and the value of 
pressure is obtained as the average of the nodal 
pressure found in each element. An improvement 
in results can be obtained by evaluating acoustic 
pressure at Gauss points inside the elements 
rather than at nodes and plot these, as is done in 
reference [6]. However, this will not totally 
eliminate the inaccuracies in the results. As seen 
in equation (10) pressure calculations require 
derivatives of acoustic potential because of the 


difficulties of this model. 

Figure 2a shows pressure distribution 
for a relatively simple acoustic field with non- 
dimensional frequency of 25. This mesh has a 
nodal density of about 11 nodes per wavelength, 
N/X, and 9021 degrees of freedom (DOF) and is 
more than adequate for finding satisfactory 
potential, p . In fact only a nodal density of 9 
nodes per wavelength is needed in certain 
regions due to the flow effect. The flow in this 
case varies from Mach .1 to .13 and the effective 
wavelength \ is given according to \ =(1"*"M)X. 
As seen, even at these very low Mach numbers 
the improper behavior of pressure is present. 
This behavior is actually so strong that the total 
number of DOF needs to be doubled before 
significant improvement can be seen. The root of 
this problem lies in the calculations of the 
potential derivatives, or the acoustic velocities, 
and the use of quadratic serendipity elements. In 
this case the acoustic propagation and the flow 
are predominately axial. This makes dpi ox the 
more significant part of in equation (10) 
when calculating pressure. Tire derivative at each 
node is calculated from the average based on 
every derivative calculated at that node, thus 
assuring continuity. But since quadratic 
serendipity elements are used dpi dx is linear in 


4 

American Institute of Aeronautics and Astronautics 


91 





x and quadratic in r while dp t cr is quadratic in 
x and linear in r. When studying the acoustic 
velocity components separately the error appears 
to lie predominately in the axial direction of 
Op l dx . The error of dpi dr in the radial 
direction, which is linear in r, does not affect the 
pressure for one main reason. dp / dr is of small 
magnitude and is multiplied with the relatively 
small flow velocity component in the r- 

direction. In addition, for this case dp / dx in the 
axial direction behaves much more smoothing 
than dpi dr in radial direction, dp! dr in the 

radial direction experiences more reflection from 
the outside wall and because of the harmonic 
appearance of the error, it is possible that there is 
some cancellation of the error because of this 
reflection since it does not appear as strong in the 
radial dp l dr as in the axial dpldx y both of 
which are linear in r and x, respectively. The 
error is of such a nature that it often appears 
more distinct in the simpler regions, not those 
heavily exposed to reflection and scattering. The 
order of the velocity components in different 
directions is at the core of the unacceptable 
behavior. A change to Lagrangian elements will 
have no effect on the order of the acoustic 
velocity terms, only on the number of terms 
involved. The increase in the appearance of error 
in certain regions, especially along the outside 
wall, occurs because the magnitude of dpidx 
increases here. The error will decrease with 
increased nodal density. 

The same acoustic field was modeled in 
figure 2b using cubic serendipity elements and 
almost the identical nodes per wavelength value. 
Again the mesh is adequate for finding 
acceptable acoustic potential. In this case there 
are only 7066 DOF due to the use of cubic 
serendipity elements rather than quadratic ones. 
There are no visible traces in the pressure 
contours of the error encountered previously. 
However, studying the acoustic velocity 
components there are small tendencies for this 
erroneous oscillating behavior, despite having 
terms of at least quadratic order. It is clear that 
the error can not be completely eliminated using 
cubic elements. This will, however, greatly 
contain it in numerous cases. The effect of the 
error on the pressure increases as the flow speed 
increases. This can be seen in figure 3a where 
Mach number ranges from .5 to .73. With the 
effective wavelength \ being more than one and 
a half of X. the value of 1 1 . 1 for N/X is far more 
than needed for satisfactory potential results. Yet 


pressure is far from acceptable in the quadratic 
case. It is important to note again that the error 
occurs in the calculation of the acoustic velocity 
and that an erroneous Vp is amplified by V^ 0 , 

not in the acoustic potential calculation (which is 
satisfactory for N l\ greater than 10). The reason 
the error appears greater as M increases is the 
product term involving V0 O in equation (9). The 
magnitude of the error in the acoustic velocity is 
not affected by flow velocity. There are similar 
tendencies for cubic elements, however not 
noticeable as seen in figure 3b. Exactly how 
much impact the presence of a moving medium 
in an acoustic field has on the modeling becomes 
clear when the flow is directed into the acoustic 
propagation. The effect of flow in the opposite 
direction makes it necessary to increase the nodal 
density to a N/X of 17 in order to account for the 
effective wavelength adjustment. The flow speed 
ranges from -.3 to -.4 and this means N/\ is 
about 10.3. Reference [6] addressed this in 
considering turbofan inlet problems and it is 
clear from figure 4a how much more sensitive 
the problem becomes. The degrees of freedom 
necessary in the quadratic case for satisfactory 
potential results have increased to 22541, yet the 
pressure contours are far from acceptable. The 
cubic model has to require 16354 DOF in order 
to have 17 nodes per wavelength and to handle 
the potential, but this is also enough to handle 
pressure. 

In cases with higher angular mode 
numbers it is generally easier to model acoustic 
pressure, as seen comparing figures 3a and 5a, in 
which the only difference is higher angular mode 
number in the latter case. The mode number 
increase appears to cause a decrease in the 
magnitude of dpidx over large regions of the 
domain, except along the outside wall. It 
basically shifts the acoustic velocity and pressure 
gradients away from the x-axis. It also causes an 
increase in dp! cr in the radial direction relative 
to dpidx . This will not increase the appearance 
of the error overall since dplci * in the radial 
direction is multiplied by the radial flow velocity 
component which is still small. 

The case involving larger contraction of 
the cross-sectional area causes the most complex 
acoustic field. This will greatly increase both the 
radial flow component and the radial acoustical 
velocity component {dpi dr in the radial 
direction increases more than that in the axial 
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Figure 3a. Pressure contours using quadratic elements for q r = 25, M; - .5, m - 1 , and N/X -11.1 


12 3 4 5 6 7 8 

0.083 0 248 0.413 0.579 0.744 0 909 1 075 1 240 



Figure 3b. Pressure contours using cubic elements for Tj r - 25, Mj - .5, m - 1, and N/X 1 1 .3 


1 2 3 4 5 6 7 8 

0.003 0.248 0.413 0.579 0.744 0.909 1 075 1 240 



Figure 4a. Pressure contours using quadratic elements for r] r = 25, Mj - -.3, m - 1 , and N/X 17 


12 3 4 5 6 7 3 

0.083 0.248 0 413 0.579 0.744 0 909 t 075 1.240 



Figure 4b. Pressure contours using cubic elements for 77 , = 25, M f - -.3, m - 1, and N/X 17 


direction), as well as reflection and scattering. 
The flow speed ranges from .1 to .24, indicating 
an N/X of 9 would be adequate for potential 
calculations. This time the magnitudes of d(p i dr 
are significant, yet because the behavior across 
the duct of d<picr is influenced by reflection 


and scatter the error might be subjected to some 
cancellations. c<p i 5x in the axial direction is 
still a problem, but now the reflection and scatter 
also make an impact in this direction for most of 
the domain. More significant, it is no longer 
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12 3 4 5 6 7 8 

0.417 0 833 1 249 1 665 2.081 2.496 2.912 3.328 



Fieure 6b. Pressure contours using cubic elements tor q r — 25, Mj - .1, m 1, and N/X 1 l.o 


totally dominating equation (10) and thus the Note that this is not why the error seems 

error it carries becomes diluted. The quadratic unnoticeable. In a smooth field larger acoustic 

model is capable of modeling this case and there velocity values will lead to the larger errors. If 

are not any significant differences when the same contour values as in previous figures 

comparing the quadratic and cubic models were used, the error would still not be noticeable 

(figures 6a and b). The area contraction will in this case, 

cause large increases in pressure magnitudes. 
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The overall tendencies seen in figures 
2a to 6b stay the same as the frequency is 
increased. Figure 7 gives an overview ol the 
necessary nodal density for acceptable pressure 
results for both the quasi-one dimensional and 
cylindrical cases. These results model the 
acoustical field in a duct with a .9 area 
contraction ratio, inlet Mach number ot .1, and 
an angular mode of 1 (, as in figures 2a and b). 
The tolerance level for what can be viewed 
acceptable has been more stringent in these cases 
for the sake of consistency. The computational 
time for two-dimensional quadratic and cubic 
elements beyond non-dimensional frequencies of 
50 and 100, respectively, becomes excessive. 



Figure 7. Nodal density vs. non-dimensional 
frequency 


The most important result is that with 
elements of cubic order the error can be 
contained enough to the point that the accuracy 
of the potential solution becomes the deciding 
factor in determining nodal density, not the post- 
processed calculation of pressure. Note that the 
benefit of going to even higher order elements, 
as is done in the one-dimensional case, is 
questionable. Very little seems to be gained. 
While overall the cubic serendipity elements 
perform better than quadratic serendipity 
elements in modeling pressure there are also 
some possible problems in its performance. 
Although nodal density might be the same in 
both cubic and quadratic models, the cubic 
elements cover a larger domain and therefore 
there is a larger gradient within this element, 
causing derivatives of greater magnitude to be 
calculated. A decrease in element size can 


correct the problem, however this will erase the 
gains made in reducing dimensionality. 

The use of cubic elements, although not 
totally eliminating the errors occurring in the 
post-processing, has contained those 
inadequacies enough to where certain high 
frequency models that previously were 
impractical to study due to dimensionality 
problems, now can be calculated in reasonable 
time. The first set of problems involves the 
acoustic field in a duct at a non-dimensional 
frequency of 50 and the impact of the area 
contraction ratio changing from .9 to .5. The 
flow velocity ranges from Mach number of .1 to 
.13 in a, from .1 to .24 in b, and from .1 to .5 in 
c, with an angular mode of 1. This model is seen 
in figures 8a to c and requires 2763 1 DOF. 

To model an acoustic field with a non- 
dimensional frequency of 100, angular modes 
ranging from 1 to 20, and the flow velocity 
varying from Mach .3 to .4 in figures 9a and b 
and Mach .25 to .9 in figures 10a and b, a 
domain of 91906 DOF is needed. Figures 9a and 
b compares angular modes. The same for figures 
10a and b, but with larger area contractions, thus 
higher flow velocities. Figures 9a and 10a shows 
the effect of changing the cross-sectional area. 
With the same degrees of freedom and a flow 
speed ranging from Mach .3 to .4, but in the 
opposite direction, the model is limited to 
satisfactory results for non-dimensional 
frequencies up to 70. The computer time for this 
model using an HP Visualize C200 with a 
specfloat FP95 of 21.4 was about 30 to 40 
minutes for the incompressible flow code and 
about 100 minutes for the acoustic propagation 
code. Each case in figures 8a to c took only 
about 15 minutes for both flow and propagation 
calculations combined. 

CONCLUSION 

It is questionable if cubic elements will 
be suitable for models having a non-dimensional 
frequency much beyond 100. It is apparent they 
can give acceptable results for nodal densities 
not much greater than those required for 
potential calculations. However, the tendencies 
of the oscillating errors are still present and this 
probably needs to be resolved before an attempt 
can be made at going below the 10 nodes per 
wavelength required for satisfactory results in 
the potential formulation. Since this error seems 
to be a function of both nodal density and the 
order of the shape functions it calls for further 
study into the use of various shape functions and 
other types of elements. 
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Figure 8c. Pressure contours for rj r = 5 0, M* — .1, m — 1, andN/X- 11*J 



123*5073 
0212 0.424 0.636 0.643 1.060 1 272 1 464 1 696 



Figure 9a. Pressure contours for r/ r — 100, Mj — .3, m — 1, and N/A 10.4 
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Figure 9b. Pressure contours for r} r - 100, Mj - .3, m - 20, and N/X 10.4 



Figure 10a. Pressure contours for q r - 100, Mj — .25, m— 1, and N/X 10.4 



1 2 3 4 5 6 7 3 

0.369 0.739 1.108 1.477 1.846 2.216 2.585 2.954 



Figure 10b. Pressure contours for rj r — 100, Mj — .25, m — 10, and N/X 10.4 
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ABSTRACT 

The boundary condition at an impedance wall in a duct with a steady mean flow requiring 
the specification of the normal component of acoustic particle velocity is examined. It is found 
that when implemented in the weak formulation of the finite element method it can be 
considerably simplified. The boundary condition would appear to require data which includes the 
tangential derivative of the tangential mean flow velocity, the normal derivative of the normal 
component of mean flow velocity, and the derivatives of the mean flow density and the boundary 
admittance along the boundary. It is shown that with suitable rearrangement the normal and 
tangential velocity derivatives can be eliminated, as can the derivatives of the mean flow density 
and admittance. The boundary condition becomes only slightly more complicated than the 
corresponding boundary condition when mean flow is absent, and is no more difficult to 
implement, requiring only local values of tangential mean flow velocity, density, and admittance 
which are already required as data for the weak formulation of the field equation. 

INTRODUCTION 

Figure 1 shows the geometry of a typical non-uniform duct section. The duct is of non- 
uniform cross section with walls S w which in general include an acoustically absorbing section 
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imbedded m an otherwise scousticslly rigid wall. Absorption characteristics of the boundary are 
given in terms of the admittance A for a locally reacting liner. The duct in Figure 1 is depicted 
as axially svmmetric, however the results obtained here do not depend on such an idealization. 
The duct geometry includes the definition of a unit normal n directed out of the fluid region, 
and therefore into the duct wall. The notional displacement of the duct wall, normal to the wall, 
is given by £ which is a function of location on the wall. In harmonic motion with time 
dependence e ' V the admittance relates this displacement to the acoustic pressure according 
to 


= Ap 


( 1 ) 


When this admittance relation is applied to acoustic propagation in ducts with steady 
mean flow it produces what appears to be a very difficult boundary condition at the admittance 
wall. Myers [1] derived the correct boundary condition which relates the normal component of 
acoustic particle velocity to the particle displacement in non-viscous flow for harmonic acoustic 
perturbations at frequency T) r as 

v-n = iti r C + F-V£ -£»•(«■ V) F (2) 


Here v • n is the normal component of acoustic particle velocity at the wall and F_ is the 
tangential mean flow velocity at the wall. Propagation in non-uniform ducts is normally modeled 
under the assumption that the mean flow and acoustic perturbation are defined by a steady flow 
potential such that V = V4> r , and by an acoustic potential such that v = Vcf) and an acoustic 
momentum equation 


p = - P r [itl r <j> + F-V(j>] 
2 


( 3 ) 


100 



A combination of equations (1), (2), and (3) produces a single boundary condition in terms of 
acoustic potential which is difficult to implement in numerical schemes. It is found that data 
required to model the boundary condition includes the derivative of the impedance and mean 
flow density along the boundary. In addition, and much more of a problem, is the requirement 
for the tangential derivative of the tangential component of mean flow and the normal derivative 
of the normal component of mean flow. These present particular difficulties because the mean 
flow data in finite element propagation models is generally obtained from a potential formulation 
for the steady flow and the required derivatives of velocity require second derivatives of the 
potential. 

The boundary condition described by equations (1), (2), and (3) has been implemented 
in a finite element scheme by Eversman and Okunbor [2], They used an approximation, argued 
to be adequate for ducts with changes in cross section which are relatively small, which ignores 
the term requiring the normal derivative of the normal component of mean flow velocity, and 
additionally ignores the effect of duct wall curvature on the calculation of the rates of change of 
quantities along the wall. The approximation for the tangential derivative of the tangential 
component of mean flow velocity is retained. The computation of this derivative is not 
considered to be very accurate. None of these approximations are thought to be significant for 
attenuation calculations in the geometries considered. 

Rienstra [3] has approached the modeling of acoustic propagation in ducts with slowly 
varying cross section by a perturbation scheme, and the analysis procedure requires the full 
modeling of the boundary condition. However, because his procedure is analytic the 
implementation of the boundary condition presents no difficulty and the issues which arise in a 
numerical model are not present. 

The motivation for the present investigation is the requirement to verify a reciprocity 
relationship which exists for acoustic propagation in non-uniform ducts with mean flow and 
absorbing linings. In order to show reciprocity no approximation in the boundary condition is 
permissible. Numerical experiments conducted with the approximate model of the boundary 
condition described in [2] suggest that reciprocity is nearly satisfied, but one is not fully 
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convinced that the small discrepancies are in the approximate model or in the reciprocity 
principle. The following work derives a model of the boundary condition for the FEM 
formulation which is exact within the FEM formalism, is easy to implement, and will replace 
the approximation in [2], Work to be subsequently reported will show that the new boundary 
condition results in a numerical substantiation of the reciprocity principle [4,5]. 

FINITE ELEMENT FORMULATION FOR DUCT PROPAGATION 

Application of finite element modeling to acoustic propagation in nonuniform ducts with 
steady mean potential flow has been previously reported [2 ]. A formulation in terms of acoustic 
potential is used to reduce the field equations to a single scalar variable. In this investigation the 
geometry of the duct and steady flow field is axially symmetric. The acoustic field is not axially 
symmetric but is represented as azimuthally periodic in a cylindrical coordinate system with x 
being the axis of symmetry, r the cylindrical radius in a circular cross section at x =0,and 0 
the angular coordinate. Solutions are sought in angular harmonics of a Fourier Series in 0 
enumerated by the angular mode number m . This reduces the solution domain to a two 
dimensional x,r plane, shown in Figure 1. The duct shape in a 0 = constant plane shows the 
surface S which defines the duct shape and could include an inner surface for an annular duct. 
Part of S includes S , which is a locally reacting acoustic treatment. 

The acoustic field is assumed to be harmonic in time at non-dimensional frequency r^. 
Geometry is non-dimensional based on a reference length generally chosen as the radius of the 
inlet at the source plane, R . Acoustic and steady flow variables are non-dimensional based on 
reference values of the speed of sound and density of the medium, , c^, generally defined at 

the plane of the acoustic source. The non-dimensional frequency is i\ r = co R / , with cothe 

harmonic source frequency. 

Reference [2] discusses in detail the finite element modeling of acoustic propagation in 
and near ducts carrying mean flow. The field equations for continuity and momentum and the 
isentropic equation of state are used in a weighted residual statement to obtain an integral 



formulation which is then written in discrete form using standard FEM procedures. In terms of 
acoustic potential the weak formulation is 



where the local non-dimensional steady flow velocity is F = V(J)^ , with (j> r the non- 
dimensional steady flow velocity potential. The local non-dimensional density and speed of 
sound are p c . The surface integral on the right hand side introduces the noise source and 
termination conditions on S Q or S L and a possible impedance boundary condition on S inside 
the duct. In the present investigation it is the impedance boundary condition which is of interest 
on S' , a portion of S. In equation (4), the weighted residuals statement, W represents an 
arbitrary weighting function selected from the class of continuous functions. In this weak 
formulation the approximation to the solution $ is also chosen from the class of continuous 
functions 

At a duct wall the mean flow is tangential to the wall and K. • n = 0 causing the 
boundary integral (the contribution to the right hand side of equation (1) related to the impedance 
condition) to become 


I b = f f p r IFV4 )-ndS (5) 

s w 

With the Myers boundary condition [ 1 ] and with V <j)^ n = 0 on the duct wall surface S w , the 
integral of equation (5) on S w becomes 

5 


103 


( 6 ) 


l b = fj{ p r FT[/r| r C - V r ‘VC - C«’(«-V)F] 

w 

The following vector identities (suggested in a similar context by Moehring [6]) are introduced 
(with account taken of the special circumstances of the present problem): 

P w V -VC = p V -VFFC - p (V-VW 

r r r r rr r r 

V*p V = 0 

r r r 

p V ■'VWC, = V-p ( 7 ) 

~r r r r 

n-V = 0 

r 

pW(n-(n-V)V r = «-(«-V)p r FTCF 

nVx(«x p r WCV r ) = V-p r PFCF - /r(irV)p r fFCF. 

With the use of the identities of equations (7), equation (6) can then be reformulated as 


I b = ff { p r C [i\W - V r -VW]}dS + f f {n-Vx(n*Wp r CV r )}dS 


( 6 ) 


Following the development of Moehring [6], the last integral can be written as a line integral on 
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the boundary T of the surface by using Stokes’ Theorem. The boundary curve r should 
enclose the portion of on which there is a non-zero admittance, but should be located where 
the admittance vanishes, as shown in Figure 2. T consists of closed curves and T 2 
circumscribed on the duct wall at either end of the duct, chosen to be outside the region in which 
the lining of finite length has non-zero admittance, that is, in the regions in which the duct wall 
is rigid. There is of course a portion of T which runs along the duct wall between T { and 
r 2 to complete the closed curve of Stokes’ Theorem, but. this curve is traversed twice, once 
in each direction, and has no net contribution. To make use of Stokes Theorem it is required that 
the acoustic field and the wall displacement be continuous on S w . Hence, if the acoustic 

treatment is of limited length imbedded in an otherwise rigid wall duct, the transition from rigid 
wall to admittance wall, as well as the variation of admittance along the treated wall, must be 
continuous. If this condition is met, Stokes’ Theorem can be cast in the form 

I f {n-Vx(nxpWCV r )}dS = f (n^pW CVJ-d? + J («* pW^V^-dT (7) 

s, r . 

The integral on the surface S ^ vanishes if the line integrals vanish. On a hard wall the line 
integrals vanish because the boundary displacement vanishes. This means that if the condition 
for the use of Stokes’ Theorem is met, then the integral of equation (6) is 


I b = f /{p r C[i\^- V-VW] }dS 


( 8 ) 


At a wall of admittance A equations ( 1 ) and (2) are used to replace the wall displacement ( 
and the pressure with velocity potential (f). The result is the new weighted residual boundary 
integral on the duct surface S w 
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j b = - f fAp 2 r {ir]W<b + FFF -Vcj) - <j>F-VFF - {V; -VW)(V r -V$)}dS (9) 

The weighted residual form of the boundary condition of equation (9) is a considerable 
simplification of the boundary condition which would result by a direct use of equations (1), (2) 
and (3). In the latter case it would be found that the derivatives of the admittance A and the 
steady flow density p are required. In addition the tangential derivative of the component of 
mean flow velocity tangential to the wall and the normal derivative of the normal component of 
mean flow velocity at the wall are required. Current implementations of the FEM formulation 
from which the steady potential flow field is obtained are not well suited for the accurate 
determination of these second derivatives of velocity potential. The modified version of the 
boundary condition neither requires data which is not directly determined from the potential flow 
model nor requires any operations which are not required in the discretization of the field 
equations (left hand side of equation (4)). 

The restriction that the admittance is continuous on the duct wall is related to modeling 
difficulties addressed by other authors. Moehring [6 ] has noted that in the acoustic potential 
formulation for discontinuous admittance variation there is no clear condition to be imposed on 
the acoustic field or wall displacement at the discontinuity. Rebel and Ronneberger [7 ] have 
shown that the condition of admittance discontinuity and the assumption of potential flow at the 
wall (no boundary layer) causes a problem with the underlying physics of the flow related to the 
absence of shear stresses. In this analysis these issues have been eliminated by requiring that the 
admittance vary continuously. In practical terms, this is accomplished by making 
“discontinuities” rapid, but continuous, variations (easily done by an appropriate definition of 
the local admittance) . One suspects that numerically this may be a non-issue, because in the weak 
FEM formulation the role of discontinuities is reduced. 
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AN ALTERNATE APPROACH 


An alternate approach to the simplified boundary condition is available which produces 
a boundary condition useful for numerical models which are not based on the weighted residuals 
formulation. The Myers boundary condition of equation (4) can be written 


p r v-fi = ir\ r C + P r 1^, ‘ V C P r ( .«•(«' V) 


( 10 ) 


The steady flow continuity equation 


V*p V = 0 

r r r 


( 11 ) 


is used to establish that 


p/ r -vc =V-p r CF 


( 12 ) 


It can also be shown that since on the duct wall n • F - 0 , 


p r C n-(n-V)V r = «-(f?-V)p r CF 


(13) 


With these results it can be shown that 


n • (fi • V) p C V = — (p C V r ) 

K r r on » 


(14) 


and 


v-p if =T (P C F).i-(pyr ) 

r r ox r r < on r 


(15) 
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Directions tangential and normal to the duct wall at the wall surface are denoted by x ,n . 

V V are the tangential and normal components of the steady flow velocity. At the duct wall 

’ r n 

V vanishes Therefore, the boundary condition on the duct wall is 

r n 


p r v-n = ix\ r p r C 


dx 


(p ,<t,> 


( 16 ) 


This form of the boundary condition, not in weighted residual form, could be used, for example, 
in a finite difference formulation. With ( replaced by equation (1) and p replaced by equation 
(3), it is found that derivatives along the wall of mean flow density, wall admittance, and mean 
flow velocity are required, however the normal derivative of the normal flow velocity component 
is not required. Equation (16) can be used in the weighted residual formulation to reproduce 
equation (9), with the same restrictions. 


CONCLUSION 

The Myers acoustic boundary condition at an admittance wall in a non-uniform duct 
carrying potential mean flow [1] has been restructured using identities of vector calculus to 
obtain a form well suited for finite element predictions of propagation. If applied without 
simplification the boundary condition would require data on the spatial derivative along the wall 
of mean flow density, the tangential spatial derivative of the tangential mean flow velocity at the 
wall, the normal spatial derivative of the normal mean flow velocity at the wall, and the spatial 
derivative along the wall of the admittance. After simplification only local values of density, 
tangential flow velocity and admittance are required. The normal component of mean flow 
velocity is eliminated completely. Implementation of the boundary condition is easily 
accomplished in finite element models. 

An alternate approach has been used to simplify the Myers boundary condition in a form 
useful for numerical modeling not based on the weighted residuals approach of finite element 
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analysis. The normal derivative of the normal mean flow velocity component at the wall is 
eliminated, however derivatives along the wall of mean flow density and velocity and wall 
admittance are retained. 

The net effect of the boundary condition on prediction of attenuation in ducts in FEM 
models has been found to be minor when compared to a former approximation introduced for 
computational efficiency (the new exact formulation is found to be even more computationally 
simple). For calculations made to validate acoustic reciprocity, the exact form of the boundary 
condition introduced here is essential, and it is found that predicted reciprocity relationships are 
accurately verified [4,5]. 
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Figure 1. An x , r section of a non-uniform duct showing essential modeling features, 
including acoustic treatment and steady mean flow. 
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in the application of Stokes Theorem. 
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ABSTRACT 

A reverse flow theorem for acoustic propagation in compressible potential flow has been 
obtained directly from the field equations without recourse to energy conservation arguments. A 
reciprocity theorem for the scattering matrix for propagation of acoustic modes in a duct with 
either acoustically rigid walls or acoustically absorbing walls follows. It is found that for a source 
at a specific end of the duct, suitably scaled reflection matrices in direct and reverse flow have a 
reciprocal relationship. Scaled transmission matrices obtained for direct flow and reversed flow 
with simultaneous switching of source location from one end to the other also have a reciprocal 
relationship. Numerical verification of the reciprocal relationships is given in a companion 
paper. 

INTRODUCTION 

The general principle of acoustic reciprocity in a medium at rest is well known and is 
derived in [1] by direct manipulation of the field equations in the case of harmonic time 
dependence. By essentially using this starting point, Eversman [2] has demonstrated and 
numerically verified reciprocal properties of the scattering matrix for acoustic modes incident, 
reflected, and transmitted in a non-uniform duct in the absence of mean flow. Moehring [3] , by 
using an approach based on energy conservation, has arrived at the same result when differences 
in definitions of the normalization of acoustic modes are considered. Moehring’s approach 
depends on a suitable definition of acoustic energy density and acoustic energy flux, which are 
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well known in the case of propagation in a stationary medium. Since the result of Moehring [3] 
depends on the normal derivative of the acoustic energy flux vanishing on the walls of the duct, it 
would appear to exclude dissipative walls. However, the classical result [1] concludes that 
reciprocal relations still hold provided the duct walls have a locally reacting impedance model 
(whether dissipative or not) where for harmonic disturbances acoustic particle velocity is 
proportional to acoustic pressure. Reciprocity based on energy conservation is more restrictive 
than necessary. 

An appropriate definition of acoustic energy and energy flux also exists for propagation 
in a compressible potential flow [4], On this basis, Moehring [3] extended his observations on 
properties of the scattering matrix to include non-uniform ducts with rigid walls and potential 
mean flow with the result that the reciprocal properties also depend on flow reversal. Godin [5] 
has studied extensively issues of acoustic energy, acoustic reciprocity, and flow reversal 
theorems in a highly generalized sense directly from the field equations. He has not specifically 
addressed the simpler case of propagation in non-uniform ducts with compressible irrotational 
flow, citing Moehring [3] as having demonstrated reverse flow reciprocity in this case [6]. 
Godin’s citations to the literature can be consulted for an extensive survey of the field. 

Here the goal is to approach the acoustic reciprocity problem in a compressible potential 
flow in non-uniform ducts directly from the field equations in much the same way as the classical 
formulation in the case of a stationary medium [1]. Furthermore, it is intended to show that 
reciprocity holds for a finite length dissipative lining imbedded in an otherwise rigid wall. A 
foundation for such a formulation in the case of uniform flows was given by Flax [7,8] in 
connection with unsteady lifting surface theory. The application to potential flows in non- 
uniform ducts given here yields a reciprocity relationship (perhaps more appropriately referred to 
as a flow reversal theorem [5]) which is in terms of acoustic potential and acoustic density 
perturbations on an irrotational compressible mean flow. It also can be given entirely in terms of 
acoustic potential perturbations. The reverse flow reciprocity formulation which is obtained does 
not begin with an energy conservation law and leads to a form similar to Moehring [3]. This 
investigation is not restricted to a duct with rigid walls. The reverse flow reciprocity theorem is 
then used to establish reciprocal properties of the scattering matrix for propagation of acoustic 
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modes in a non- uniform duct with acoustically absorbing walls. 

In the present paper the theoretical frame work is established for the reciprocity (flow 
reversal) relationship, and it is specialized for examining the reciprocal relations which exist 
between duct modes propagating in compressible mean flow. In a companion paper [9], the 
results derived here are substantiated by numerical experiments based on a finite element 
simulation using a new implementation of the boundary condition introduced by the presence of 
acoustic treatment, the subject of a second companion paper [10]. 

ACOUSTIC PROPAGATION IN A COMPRESSIBLE POTENTIAL FLOW 

An extensive discussion of both linear and non-linear formulations for acoustic 
propagation in potential flows has been given by Campos [11]. His work includes a number of 
citations to previous work and is the basis for contributions directed mainly toward analytic or 
semi-analytic solutions for propagation in ducts (see for example [12,13]). The investigation 
reported here has been part of the development of numerical modeling methods for acoustic 
propagation in non-uniform ducts and therefore the final form of the governing equations is 
specialized for that purpose. 

The acoustic field equations are obtained by the consideration of unsteady perturbations 
on a steady compressible potential flow. Accuracy in calculation of both the steady and unsteady 
flow fields is necessary for computational verification of the theoretical results obtained. The 
starting point for the formulation of both the steady mean flow and the acoustic perturbation 
consists of the mass and momentum equations and the energy equation in the form of the 
isentropic equation of state: 


+ V-(pF) - 0 C 1 ) 

dt ^ 


— + (F-V) V = ( 2 ) 

dt p 
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( 3 ) 


p , p , F are fluid properties pressure, density, and velocity, at this point in dimensional form. 
p Q , p Q are reference values of pressure and density. It is assumed that the mean flow and 
acoustic perturbations are irrotational and that a potential <f) exists such that V = V <p . Acoustic 
perturbations are assumed on the steady mean flow such that 4> = <|> r + <l>,P=P r + P and 
p = p r + p. The linearized continuity equation is 

^ + V-(p r V4> + pV<j> r ) = 0 (4) 


The linearized momentum equation, for irrotational acoustic perturbations, is 


p = -fl ( + V<i> • v<j) ) (5) 

H c 2 at r 

r 

This is used to replace p in equation (4) and the linearized equation of state, 

p-c r 2 p (6) 

is used to produce an alternative form of the momentum equation in terms of acoustic pressure, 

p = -P r (^ + V<|> r -V<j>) (7) 

Equation (7) is used to post-process the field solution for (J) to obtain the acoustic pressure field. 
The acoustic particle velocity and acoustic velocity potential are related according to 


v = V<j) 


( 8 ) 
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The perturbation process also produces the conservation equation for the steady flow 

V'(p r V(J) r ) = 0 (9) 

and the steady flow momentum equation in terms of the speed of sound 

c r 2 = \- ^ ^[V4yV(j) r -M 2 ] (10) 

and in terms of the steady flow density 

l 

P,*[l (U) 

Equations (4) through (1 1) are now in non-dimensional form where $ is the acoustic potential, 
<j> r is the local mean flow (reference) potential, p is the acoustic density, p r is the local mean 
flow density, and c r is the local speed of sound in the mean flow. All quantities are made non- 
dimensional by using the density p_ and the speed of sound at some point, in this case the 
radius at the plane x = 0 . Stagnation conditions could also serve as the reference. A reference 
length R is defined as some characteristic dimension at the plane x = 0 . In the case of a 
circular duct the reference length is the duct radius at x = 0 . The acoustic potential is non- 
dimensional with respect to c^R, and the acoustic pressure with respect to pj 2 . Lengths are 

made non-dimensional with respect to R. Time is scaled with . In the case of harmonic 

time dependence this leads to the definition of non-dimensional frequency v\ r = (x>R/c^. wis the 
dimensional source frequency. Af is the Mach number at the reference point. 

Equation (9) is the field equation for the calculation of the compressible potential mean 
flow. Equations (10) and (1 1) are subsidiary relations that can be used in an iterative solution 
which at each stage uses a density field derived from the previous iteration step. VcJ)^ , c ^ , p^ are 
required data for the formulation of the acoustic problem. 
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REVERSE FLOW RECIPROCITY PRINCIPLE 

The application of a reciprocity relationship for acoustic propagation in potential flows in 
non-uniform ducts has lagged behind the exploitation of the comparable results for propagation 
in a quiescent medium, not because of the difficulty in posing the principle, but probably because 
of the difficulty in producing solutions which could be used to test it. Numerical solutions for 
duct propagation using the Finite Element Method (FEM) are now achievable [14] and provide 
the capability of determining the complete acoustic field in a duct (and in the far field of a duct 
of finite length) as well as the scattering matrix for a non-uniform duct inserted in an otherwise 
uniform duct of infinite length. This provides opportunity for testing the reciprocity principle and 
suggests the development of such a principle for benchmarking of FEM calculations. 

The intent here is to approach the reciprocity principle independent of considerations of 
energy conservation. A counterpart exists in the literature of unsteady lifting surface theory in the 
Reverse Flow Theorem of Flax [7,8] . A reciprocity principle for acoustic propagation in non- 
uniform potential flows in ducts can be obtained by an extension of the formulation of Flax. 

Consider the volume Q shown in Figure 1, which is the interior of a nonuniform duct of 
arbitrary cross section. In examples the duct will be assumed axisymmetric (circular or annular), 
but the principle derived is independent of the duct cross section. The duct walls can be rigid or 
locally reacting. The unit normal n is directed out of the volume at each surface. The source 
plane S is where the acoustic source is specified and the exit plane S g terminates the duct and 
may have a reflection matrix specified. For computations the exit plane will be assumed non- 
reflecting. A typical computational problem would seek to specify the acoustic field within the 
duct and the scattering matrix at the source plane for incident acoustic modes. Equations (4) and 
(5) specify the acoustic field within Q subject to appropriate boundary conditions on S, the 
surface of Q . 

Let <j) e lT]r ‘ be a harmonic solution for the acoustic velocity potential for the case of a 
mean flow specified everwhere in the duct by its reference Mach number M = V 4> r and with 
specified boundary conditions. Let c[> 2 e be a second harmonic solution for exactly the same 
duct with different source conditions, but with the flow reversed, -M f - - Vcj)^. It is important 
to note that in reversed flow the reference density p r and reference speed of sound c r are 
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unaltered. Because of equation (4), in the case of a harmonic source at non-dimensional 
frequency r\ r , it follows that 


f ff(4> 2 [ir] r p l + V-(pV<j) l + V^p,)] - <t> 1 [fri r p 2 + V'(p r V(j) 2 - V(f) r p 2 )]}^Q = 0 (12) 

Q 

With application of the divergence theorem equation (12) is reconfigured to 

//{<i> 2 [(p r v<iy V(j) rPl )] - 4>,[(p r v<j> 2 - V(|) r p 2 )]}-«^ 

5 

- Iff + Vcj) 1 -V(j) r p 2 - /T lr (* lPl - 4)^,)}^ = 0 (13) 

0 

The acoustic density in the two solutions is defined according to 


-- sO'Vfrj + V^-Vcj),) 


-£0\4> 2 - Vcj) r -V4) 2 ) 


Equations (14) are used to eliminate the remaining volume integral in equation (12), leaving the 
reciprocity principle in terms of acoustic density and potential in a form convenient for 
subsequent development 


//{^[(P.Vct), »V4> r p l )l-<|. 1 [(p,V* 2 -V* r p 2 )]>- B - fl B = 0 (15a) 

An alternate form, in terms of the acoustic pressure, obtained by using equation (7), is 

I J{(_! - V<J> r -V4> 2 ) Kp r V<|> 1 + V^Pj)] - (— + V4) r -V(l) 1 )[(p r V(j) 2 - V<b r p 2 )]}-ndS = 0 
s p r r r P ' (15b) 


In establishing reciprocal relationships for the scattering matrix it is important in this 
development that Equation (15a) or (15b) have contributions only on the source and exit planes, 
S and S . This requires that the integrand vanish on the duct walls. For a duct with rigid walls 

s e 

7 


119 


this occurs because Vtjyn = 0 , V<ty» = 0 , and VcJyH = 0 on the walls. In the case 
when mean flow is absent , V 4^ vanishes and is constant, leading to 


ff(P 2 v } ~ P^ij-ndS = 0 ( 16 ) 

For a locally reacting wall admittance in the absence of mean flow, acoustic pressure is 
proportional to particle velocity. This makes the integrand vanish on the walls of the duct, and 
equation (16) has contributions only on the source and exit planes. Reciprocal properties of the 
scattering matrix are therefore valid in the absence of flow for normally reacting impedance 
walls. This is in spite of the fact that acoustic power is not conserved. 

When mean flow is present and the walls are normally reacting further examination is 
required to establish that equation 15(a) has no contribution on the duct walls. Myers [15] has 
shown that the boundary condition at a normally reacting acoustic lining which relates the 
boundary displacement ( to the component of the acoustic particle velocity normal to the 
undisplaced surface is 

V(p-n = /n r C +M/VC -C«-(«'V)M (17) 

With V(b ’H = 0 on the duct wall surfaces 5 , lined or unlined, and M = Vtj)^ the integral 

of equation (15a) on S ^ becomes 


I w = / /{pA[”Vn - C,«-(«-v)M] 

K 

- [iri r C 2 -M-V C 2 + C 2 n-(n-V)M r ]}dS 


( 18 ) 


The following vector identities are introduced (with account taken of the special circumstances of 
the present problem): 


p r 4 )M-vc = p r A?-v<K - p r C^/V4> 
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V-p M r = 0 


pM r -V4)C = v-p r 4 >CM 

«-M = 0 

r 

p^4 )(>r(/rV)M = «-(«-V)p r <j>£M r 
n-Vx(«xp4)CM) = V- P r 4>C M r - «-(«-V)p r (|)CM 
Equation (18) can then be reformulated as 

/„ = //{ Pr c,[inA -w-v<t> 2 ] -p/jftiA *H r -^ I ])ds 

s „ 

+ If {^' Vx (^ x P r ^2^A^ + n-Yx(nxp r (P i ( 2 M)}dS 

S w 

(19) 

The case considered here is an admittance wall imbedded in an otherwise rigid wall. Therefore, 
as shown in figure (2), acoustic treatment extends less than the full length of the duct. This is 
physically realistic, and is also consistent with the type of reciprocity relation between duct 
modal amplitudes which is sought. Following the development ofMoehring [16], the last 
integral can be written as a line integral on the boundary V of the surface by using Stokes’ 
Theorem. Since the surface consists of the wall of the duct (in general of varying cross 

section), r Is chosen to consist of closed curves and T 2 circumscribed on the duct wall 
outside of the region of the lining, and therfore where the duct wall is rigid. There is also a 
portion of V which runs along the duct wall between and T 2 to complete the 

closed curve of Stokes’ Theorem, but this curve is traversed twice, once in each direction, and 
has no net contribution. To make use of Stokes’ Theorem it is required that the acoustic field and 
the wall displacement be continuous on S w . Hence, since the acoustic treatment is ot limited 
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length imbedded in an otherwise rigid wall duct, the transition from rigid wall to admittance 
wall, as well as the variation of admittance along the treated wall, must be continuous. If this 
condition is met, Stokes’ Theorem can be cast in the form 

[f{n-Vx(nxp<PZM)}dS = f (fixp^C M)-dT + f (nx p^C M)-dP (20) 

l r . 

The integral on the surface S w vanishes if the line integrals vanish. On a hard wall the line 
integrals vanish because the boundary displacement vanishes. This means that if the condition for 
the use of Stokes’ Theorem is met, then the integral of equation (19) is 


/. = //{p,C 1 [iVt’ 2 -W-V<i) 2 ] - P,C 2 [nvt>, + Ayv^]} dS (20) 

w 

With equation (7) this becomes 

W/< (21) 

s. 

At a wall of admittance A there is a relation between pressure and wall velocity which is 
frequency dependent and of the form 

in r C -A P < 22 > 

The integral I vanishes and equation (14) has contributions only on the portion of the surface 
area S which corresponds to duct cross sections beyond the impedance wall, on the surfaces 
S And S . With appropriate restrictions on the impedance wall, the reciprocity principle is 

s e 

therefore unchanged for hard and soft wall ducts. 

The restriction on the impedance wall is interesting. If admittance is indeed discontinuous 
along the wall, thenC has to be discontinuous. Moehring [16 ] has noted that for discontinuous 
admittance variation there is no clear condition to be imposed on the acoustic field or wall 
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displacement at the discontinuity. Rebel and Ronneberger [17 ] has shown that the condition of 
admittance discontinuity and the assumption of potential flow at the wall (no boundary layer) 
causes a problem with the underlying physics of the flow. What is known is that a numerical 
procedure such as the FEM restricts the solution for acoustic potential to continuous functions, 
and the lining displacement never appears in the final field equations and boundary condition. 
Thus one suspects that the restriction to continuously varying admittance is not a critical issue in 
numerical comparisons, but it may be in comparisons with experiment. 

Equation (14) can also be written entirely in terms of acoustic potential 


ff<t> 2 {p r ( v <t>i - -^A' vc M - 

S C r r 


S °r r 


(23) 


where the area S is now understood to include only the source and exit planes and 
S . 


APPLICATION TO A NON-UNIFORM DUCT 

The discussion here is presented for a duct with a straight x axis, but is more 
complicated only in the notation if the axis is not straight. Figure 2 shows a representative non- 
uniform duct with cross section S(x) defined on 0 <x<L , S(0) = S Q , S(L ) = S L There is a 
steady mean potential flow in the duct M(x,r), defined as 




1 — Vd> 

c r {x,r) 


(24) 


A/ is a non-dimensional velocity based on the reference speed of sound c ^ . A/(x,F)is the 

r 

local Mach number defined in the usual way and c(x,r) is the local (dimensional) speed of 
sound. (x,r) denotes a point in a cross section at x in a coordinate system appropriate for the 
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duct geometry. The corresponding reversed flow is - M . The non-dimensional frequency r\ 
can also be defined locally (local speed of sound but with reference length R) according to 


r| is a local non-dimensional frequency based on the local speed of sound. In terms of local 
Mach number and non-dimensional frequency, equation (23) becomes 


|| P r 4>2 { - ir\M§ { }-fidS 

- || { [V4> 2 -M(M‘ v 4> 2 )] + ir\M<b 2 )-ndS = ° 

h* s i 

(26) 


It is assumed that at the inflow end of the duct at x = 0 the steady mean flow is uniform on the 
cross section with M{x,f) = M Q i , and at the outflow end x = L the steady mean flow is 
uniform with M(x,r) = M.i. Reference density p and p and local non-dimensional 
frequency q Q and r} L are defined similarly. The assumption of uniform conditions implies that 
the inflow and outflow planes are well removed from the non-uniform region of the duct. In 
computational examples it is found that for ducts with circular or annular cross sections uniform 
inlet and outlet ducts of length two duct radii ahead of and beyond the non-umformity are 
sufficient. At x = 0 the outward unit normal n = -i and at x - L the normal is n - i . With 
these observations the reciprocity principle of equation (26) becomes 


12 


124 


a<t>, 


a*. 


8(j), 

= p r //< MO 


(27) 


„ act>, 

-p r J/<M(i AA<M 




CIRCULAR DUCT EXAMPLE 

In the regions of uniform flow at the ends of the duct at x = 0 and x - L acoustic 
potential can be approximated by an N term eigenfunction expansion in terms of duct modes 
(Figure 2). In the case of a circular duct the expansion can be expressed at x = 0 in vector- 

_ ]m g 

matrix form for angular dependence e as 


(b (x,r,6) = $ O) 


e'(x) 


a \e ~ ,m% + 

m 


«.w]kHK)*'" e (28) 


The derivative is 

84) 


ox 


(x,r,0) = $(r) 


«:(x) 


ik* {a + \e + 

x 1 m 




«'(*) 


■ ik ‘ 


K 


e'"” 0 (29) 


[$(r) j is a 1 x /7 row matrix of duct radial modes, the same for both right and left propagating 

are N N diagonal matrices with typical elements s 


modes. 

e m 

and 

e'{x ) 

m v 

-ik* 

and 

-ik~ 

are 

X 


X 
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and \a 


are N* 1 vectors of modal amplitude coefficients for right (positive x) and left 
(negative x) modes. A similar expansion with modal amplitudes |& m j and |& m | applies at 
x = L . The axial wave numbers are given in the nominal flow for modes which are cut on by 


/ k ± ' 

1 mn 


\ ' l 




/ \ 
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-M± . 

N 

1 - ( 1 - M 2 ) 

k a 

mn 


<N 
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(30a) 


and for modes which are cut off by 
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mn 
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(30b) 


In reversed flow for cut on modes 
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and for cut off modes 
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(31b) 


The non-dimensional frequency T) is based on the local speed of sound and the reference 
radius, k o are eigenvalues determined from the uniform duct eigenproblem [18] in a uniform 
duct with local radius possibly different than the reference radius. In the case of the circular duct 
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they are determined from Jj{ k^o) = 0,with J m being the Bessel function, m is the angular 
mode number and n , 1 <. n < N, is the radial mode number, a is the ratio of the local duct 
radius to the reference radius, o=R,/R. ot\ is therefore non-dimensional frequency based on 
local speed of sound and local duct radius. At x = 0 where the reference radius and c m are 
defined, o = 1 and p = q r . At x = 0 the diagonal matrices [^(0)] and [<? m (0)J 
become identity matrices. At x = L [<(Z)j and [e m (L)\ can be absorbed into the 
amplitude coefficients. The convention on the sign choice in equation (30) corresponding to 
k * is that the positive sign is chosen if the radical is real and the minus sign is chosen if the 
radical is imaginary. The opposite choices are made for k~ . k* then corresponds to waves 
propagating in the positive x direction (except for the possibility that with Mach number negative 
some propagating waves may appear not to propagate in the positive x direction due to 
convection) and to cut-off modes decaying in the positive x direction. The opposite 
interpretation applies for k 

mn j 

With these observations, deleting the implied dependence on the mode number m , and 
taking the reversed flow solution to have angular dependence e ,m0 , it is possible to express the 
solution <j> and its counterpart reversed flow solution 4? 2 as 

i> l =[«']K) e ' , ’ e T*]KK , " e (32: 



The corresponding reversed flow solution is 



(33) 


(34) 


34> 2 

dx 
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It is important to note that the eigenvalues and eigenfunctions are independent of angular 
dependence e~ ime or e 0 . The choice of angular dependence in the two solutions eliminates 
the angular dependence in the integrals of equation (27). The physical implication is that of a 
spinning mode which has the same vector sense with respect to the flow direction. In 
calculations, scattering coefficients do not depend on the sign of m. 

In carrying out the integrals of equation (21), the notation 




j 


j. 


(36) 

(37) 


is introduced. 


J, 


and J 


are N x N diagonal matrices resulting from the orthogonality 
of the duct eigenfunctions at x = 0 and x-L. For a circular duct the eigenfunctions 
are Bessel functions of the first kind of order m . In numerical implementations it is convenient 


to venerate the Bessel functions, solve the related eigenproblem, and generate 


J 


and 


J r 


using an FEM formulation. 

With the eigenfunction expansions of equations (32)-(35) the integrals of equation (27) 


can be written 


34 >! 


P / / ' M ^~dx " iX] 0 M 0^^ dS 

= P ro [(l - M 0 2 )({<} r + {a;} T )[J 0 ]([~ik; J«> + H*' JK)) 

- iq 0 M 0 ({a;} T + {a;} T )[J 0 ]({a;} + {a;})] 




P rjf^ 1 {(1 " + iV[ 0 M 0^ dS 


(38) 


= P, o [ ( 1 - M 0 J )({a,*> r t {a 1 '} r )[^K[-«; o2 lf<} * [->'*; 

* in 0 M a ({a;} T * {a, } r )[J 0 K {a 2 *> * {Oj-})] 


(39) 
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34 ), 

p r //4> 2 {(1 - ^ dS 


= P r [(1 ~ M £ 2 )({6 2 ‘} r + {b 2 } T )[J l ]([ lk XLJ ]{b l }+[ ik x Lj ]{b j }) 

- iT] L M L ({b*} T + {b 2 } r )[^ £ ]( {b { } + {by })] 


( 40 ) 


34), 

P r //4>j{( 1 -M 2 L )^-i\M L ^}dS 


P ^[(i - M 2 L )({b;} T + {b;} T )[J L ]([-ik; L ]{b;} + [~ ik x J{ b i^ 
+ zr^M^ ({6^} r + {by } r )[^]( {b 2 } + {b 2 })] 


(41) 


The notation k* designates the axial wave number evaluated at x - 0 for the nominal flow 

direction, k* corresponds to reversed flow. kT designates the axial wave number 
x 02 LI 

evaluated at x = L for the nominal flow direction, k* corresponds to reversed flow. 

12 

In equations (38)-(41) amplitude coefficients { a ± ) are associated with an eigenfunction 
expansion at x = 0 and { b *} correspond to x = L . With the use of equations (30a) and 
(30b), introduce the following definitions for the nominal flow direction for propagating modes 


( k real) 


a\ = p r [~i{l ~ M 2 )k*^ - ir\M] = -i P r fl^ 1 - (1 “ M 1 ) 


( \ 

KO 


(42) 


a‘, =p r [-i(l -M 2 ) AT - it\M] = /p^J 1 - 0 ~M )^ — 


KO 


(43) 
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and for cut-off modes ( fc. complex) 


a\ = p r [-/(l - M )k* ^ - ir\M] = -p r ^ 
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(1 - M 2 ) 

KO 

- l 
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(44) 
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p [-i(l - M )k~ - ir\M] = p r P 

r * 1 \ 
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2 

(1 - M 2 ) 

KO 

- l 




(45) 


Analogous definitions can be introduced in reversed flow using equations (31a) and (31b). For 
example, for a propagating mode 


a 


= p r [-z'(l - M 2 )k ^ + ir\M] = -/ P r r l* 1 ~ (1 “ M 2 ) 


KO 

or l 


A 


(46) 


) 


It becomes apparent that the definitions do not change in reversed flow so the conclusion is made 

that cf = cf = a' and cf = cf = cf for both propagating and cut-off modes. 

2 1 2 1 

Furthermore, it is apparent that of = - of and these are to be evaluated at x = 0 and x = L 
as required. 

Equations (36)-(41) can be rewritten in the form 


S4>, 

pJ/<MO 

\ 

(47) 


18 


130 



3<t>, 

(48) 

= «>V 0 ][«'„]K> ♦ <a 1 ') I 'KK“‘ 0 ){a 2 '} 

* K> I V 0 H«‘o]{‘ , 2> * K>’V ( ,H%H a 2 '> 

act>. 

P //<t> 2 id 

= {kTt’VjKJW * 

* {6j'} r [J' i H«' I ]{*,'} * < 49 > 


a<i>, 

p r // ^i 10 ‘^“ar “ 

^z. 

= {*,-} r [^][cx*J{6 2 -> M4,*m t ]KJ(* 2 ‘} 

+ {b { } r [*/ £ ][« £ ]{^ 2 } + {^i } r [^,][ a J {* 2 ) 

(50) 


The diagonal matrices [a^] and [a ± £ ] are have elements defined by equations (42)-(46). 

Elements ct are evaluated at x = 0 and cc are evaluated at x = L and the distinction 
o ** 
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between nominal flow and reversed flow disappears. At x = L the modal amplitudes are b n 
and b ' . Modal amplitudes a' n and a~ and b ' n and b n are related by the acoustic 
potential scattering matrix according to 



( 51 ) 


where the scattering matrix is defined as 


[R] [ f ] 
[T] [*] 


(52) 


Contained in [5] are the usual reflection matrix [i?] and transmission matrix [T] for acoustic 
modes incident at x = 0 and reflection and transmission matrices [R] and [f] for modes, 
incident at x = L . There will be a scattering matrix ] for nominal mean flow and a second 
one [S' ] for reversed flow. The relationship between [5^] and [S^] can be obtained using the 
reciprocity theorem. 

When the integral evaluations of equations (47)-(50) are used in equation (27) there is 
considerable simplification due to the fact that cc Q = - cc Q = - a Q and a L = ~ a i = ~ a i • 
The diagonal matrices [a ] and [a £ ] are constructed by evaluating equation (42)- (45) at 
x = 0 or x = L for each acoustic mode included in the acoustic potential expansions of 
equations (32) or (35). [J Q ] [a Q ] and [JJ [aj are diagonal and therefore equal to their 
transpose. The result of the simplifications is 


K} r [/ 0 ] [«„]{<} ♦{ b;} T [J L ][a L ]{b -} = 

{a-} T [J 0 )[a 0 ]{a;} + {b;} T [J L ][a L ]{b;} 


(53) 
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Equation (53) can be written in partitioned form 
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(54) 


Equation (54) is rewritten by introducing the definition of the scattering matrix from equation 
(51) and by using the definition 


WKl 

W 


[</][>] 


(55) 


The result is 


< 


a t 


[S'/ [./][«] 




[5 2 ] r [J][a] 



(56) 


Equation (56) reveals that 

[S { f[J ][<*] = [J] [a] [S 2 ] 


(57) 


or 


[j] t a ] [£,] = ([^[«][5 2 D r 


(58) 


Equations (57) and (58) show that a weighted version of the nominal flow acoustic potential 
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scattering matrix and similarly weighted version of the reversed flow acoustic potential scattering 
matrix are transposes of one another. In terms of the acoustic potential reflection and 
transmission coefficient matrices the result is 


[*,f [/,][«„] = [•',] K 1[« 2 1 

(59) 


(60) 

[r,] r [J t ][« i ] = [4] [«„][?*] 

(61) 

[ 4 ] [43 = [41 1 4 1 

(62) 


The reciprocal relationships of equations (59)-(62) involve acoustic potential reflection and 
transmission coefficient matrices, with diagonal elements representing reflection and 
transmission coefficients in the incident modes (here referred to as direct reflection or 
transmission) and off diagonal reflection and transmission coefficients from the incident mode to 
another mode. Equations (59) and (60) show that direct acoustic potential reflection coefficients 
are invariant in reversed flow. The transmission coefficient matrix pairs [ 7^ ] , [T^] and 
[f ] , [f 2 ] are not reciprocally related but the pairs [TJ , [fj and [fj , [7^] are related 
by equations (61) and (62). These results for acoustic potential reflection and transmission 
coefficient matrices are more interesting than those obtained for acoustic pressure reflection and 
transmission coefficient matrices in the absence of flow [2] because they identify a relationship 
between reflection coefficient matrices in nominal and reversed flow which includes the 
observation that the direct reflection coefficient (in the incident mode) is invariant to flow 
direction. 
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RECIPROCITY IN TERMS OF ACOUSTIC PRESSURE 

The entire development to this point has been earned out in terms of acoustic potential 
because the field equations (4) and (5) favor this formulation. Equation (7) provides a direct 
relationship between acoustic pressure and acoustic potential which can be used to restructure the 
reciprocity results in terms of acoustic pressure modal amplitudes. In terms of local non- 
dimensional frequency and local Mach number in a uniform section of duct with uniform flow, 
equation (7) can be rewritten as 

p = -p c (zr|(j) +M~) (63) 

v r r dx 

With equation (63) a connection between acoustic potential modal amplitudes and acoustic 
pressure modal amplitudes can be established. Consider an acoustic mode propagating in the 
uniform section with the axial wave number given by equations (30 a,b) or (31 a,b). By referring, 
for example, to equation (33), the pressure amplitude can be found in terms of the potential 

amplitude in nominal flow as 
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± 
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-inp,<7 (i 


k~ 

M—) 4>j ±= 

I 



(64) 


and in reversed flow by 



-i'flP r c r (1 


k ± t 

M— )**, = 


(65) 
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These relations between acoustic pressure modal coefficients and acoustic potential modal 
coefficients are evaluated at x = 0 and x = L to produce transformations between acoustic 
potential modal amplitudes { a ± } , { b ± } and acoustic pressure modal amplitudes 
{**},{**> : 
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Po 



[S ; ] 


vl 


d i 


(66) 



(67) 


( 68 ) 


(69) 


The 2 N x 2 N diagonal matrices [B ± ] and [5 1 ] have coefficients defined by equations (64) 
and (65) for each mode, evaluated at the appropriate end of the duct, arranged along the 
diagonal. Equations (66)-(69) and equation (51) provide a relationship between the scattering 
matrices for acoustic potential and the scattering matrices for acoustic pressure. 


[S,] = [ST‘ [S,][a*l 


(70) 
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[SJ 5 [5^] are the scattering matrices in nominal and reversed flow for acoustic potential modal 
amplitudes and [Sj ] , [SJ are the scattering matrices for acoustic pressure modal amplitudes. 
Equations (66)-(71) and equation (56) are used to arrive at the reciprocity relationship for 
acoustic pressure modal amplitudes in terms of the nominal flow and reverse flow acoustic 
pressure scattering matrices: 

[£*]-' [s/tnw [s'] = [smsiMtsjnsT 1 < 72 > 


The four reciprocal relationships are: 

[P 0 t‘ [s.fKi [“oi cp»] - [P®] t- 7 .] [“») [*2] tP 0 ] 1 

[p-r 1 ki ip;i - [PH ty Ki i tp ;] 1 
[Pj ]' 1 [T;f[/ t ] ip;i - [p;i Ki t« 0 ] fr] [PT 1 

cp-r‘ [TfiJ,] [«,] [PJ] = [PP 1^1 [«j fr 2 ] [P 0 ]'' 


.(73) 

(74) 

(75) 

(76) 


The reciprocal relationships of equations (73)-(76) involve acoustic pressure reflection and 
transmission coefficient matrices, with diagonal elements representing reflection and 
transmission coefficients in the incident modes (here referred to as direct reflection or 
transmission) and off diagonal reflection and transmission coefficients from the incident mode to 
another mode. From equations (73) and (74) it is seen that in terms of acoustic pressure modal 
amplitudes direct reflection and transmission coefficients in nominal flow and reversed flow are 

not invariant but are simply related. Transmission coefficient matrix pairs [7^] , [7^] and 

[fj , [t,] are not directly related but pairs [7j , [£j and [?, ], [T 2 ] are related by 



equations (75) and (76). The result in the case of a plane wave incident is particularly simple and 
will be given in the following section. 

Because of the way in which equations (73)-(76) were developed, the weighted, or scaled, 
matrices on the right and left hand sides are equivalent to their counterparts in equations (59)- 
(62). This convenient definition of scaled pressure reflection and transmission matrices makes 
them numerically equal to their scaled acoustic potential counterparts. 


ACOUSTIC POWER CONSIDERATIONS 


Reciprocal relations between scattering coefficients in uniform duct sections bounding a 
non-unifom section obtained above are of theoretical interest and are also a useful tool for 
benchmarking numerical models of duct propagation. Also of interest in this respect are acoustic 
power transmission and reflection characteristics of the nonuniform section. For a rigid wall duct 
acoustic power conservation provides a valuable benchmark and for an acoustically treated duct 
acoustic power transmission calculations are required to assess performance. In this section 
acoustic power formulations appropriate for the propagation model are obtained. 

There are two commonly used types of acoustic intensity formulations in moving media 
[19], The Type I definition due to Morfey [4] is valid as part of a conservation law in non- 
uniform ducts for compressible potential flow, while the Type II formulation due to Ryshov and 
Shefter [20] is valid as part of a conservation law only for uniform flow, and therefore only for 
uniform ducts. In the following development the Type I definition of acoustic intensity is used in 
the uniform flow sections on either end of the nonuniformity to obtain acoustic power 
expressions. 

The Type I acoustic intensity is defined in non-dimensional form as the time average 
acoustic energy flux 



pv + 


p r c r (M-v)v + 


— Mp 2 

Pr C r 


M(M-v)p 


(77) 
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where M is the local Mach number of the flow and p r , c. are local values of the non- 
dimensional mean flow density and speed of sound in the mean flow. For propagation m uniform 
flow in the direction of the x axis equation (77) simplifies to the scalar form 




+ M 2 )pu 


+ p c Mu 2 + 

~ r r 


l 

p c 

r r r 



(78) 


Acoustic power is obtained by integration over a cross section. For the circular cross section 

= f f ( ( 1 + M 2 ) pu + p c r Mu 2 + Mp 2 | dS (79) 

A ref pyi s \ r r Pr ° r I 

A modal expansion as given by equations (32) and (33) is used to obtain an acoustic power in 
terms of acoustic potential modal amplitudes. Because of orthogonality of the duct 
eigenfunctions and because the duct eigenfunctions are the same for propagation with and against 
the flow, the evaluation of acoustic power is considerably simplified. The result can be cast in 
the form of a “power matrix” 
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ref 


P 

1 OO 00 




(80) 


where, for example, a ’ denotes the complex conjugate of a . The power matrix [P] is 
structured in diagonal blocks due to orthogonality of the acoustic eigenfunctions. The diagonal 
blocks consist of power coefficients for positive and negative acoustic modes. The off-diagonal 
blocks represent power due to interaction of positive and negative modes with the same 
eigenfunction. If the amplitude coefficients are for acoustic potential, the power coefficients are 
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There is no interaction power for cut on modes. 

For cut off modes 
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and the acoustic power transmission coefficients are 

P ** = ?'* = 0 
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P'* = — r] 2 p c J 
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(85) 


For cut off there is only power related to the interaction of positive and negative modes. It is noted 
that the acoustic power coefficients based on acoustic potential amplitude are invariant to the 
direction of the flow. 
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In the case the amplitude coefficients are for acoustic pressure, the power coefficients can 
be obtained by making use of equation (64) which relates pressure amplitudes to acoustic potential 
amplitudes. The results are 
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The conclusions drawn from this formulation in terms of acoustic pressure amplitude coefficients 
are similar to those obtained from the acoustic potential form. The acoustic power coefficients for 
the case of cut on modes are 
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For modes which are cut off, the acoustic power coefficients are 
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The power coefficients in terms of acoustic pressure amplitudes depend on the direction of the mean 
flow, but conclusions drawn about the vanishing or non-vanishing of the power coefficients for cut 
on and cut off modes are the same as in the case of acoustic potential amplitudes. 

The common result of both formulations is that elements of the diagonal blocks of the power 
matrix vanish for cut off modes and elements of the off-diagonal blocks vanish for cut on modes. 
For cut on modes there is no power contribution due to interaction of positive and negative 
propagating modes. For cut off modes there is no power contribution due to individual positive and 
negative modes but there is power due to the interaction of positive and negative modes. 

Acoustic power at x - 0 is written in terms of the power matrix and the amplitude 
coefficients (either potential or pressure formulation) 


n o = ^"> r [C 0 H«U + {a*V[C 0 HO 

+ {“'~} r [ p ™ 0 ] {a > + {a ' }T[P ™ 0 ]{a } < 89) 

With the source considered to be at x = 0 , the reflection matrix can be used to replace reflected 
modal amplitudes to yield 
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n„ = (a-’l'trH..-) * 

(90) 

This result requires interpretation. If only cut on modes are included among the incident modes a ' 
then there is no interaction power and the net power is accounted for by the sum of incident power 
and reflected power, represented by the first and last terms. If among the incident modes there are 
ones which are cut off, then there is an additional component of net power due to interaction of 
incident and reflected cutoff modes, represented by the middle two terms. In this case the net power 
is not conveniently partitioned between incident and reflected contributions. 

Acoustic power at the exit end of the duct, x - L, is similarly written in terms of modal 

amplitudes 

n L = {b"f[P~ ]{b + } 

(91) 

The termination is assumed to be reflection free, so only right modes are present. The modal 
amplitudes are related to the source modal amplitudes via the transmission matrix. This yields. 

n £ = {a* + } r [r;f [ p;; ][f,](^} 

(92) 

The form of equations (89)-(92) is the same in direct or reverse flow, however for multi-modal 
propagation no simple relation between acoustic power in nominal and reversed flow appears to 
exist, nor is there a simple relation when the source location is reversed. As will be shown m a 
companion paper, there are simple relationships for power under flow reversal and source reversal 
for a one dimensional model of propagation valid at low frequencies. 
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For a rigid wall duct energy conservation requires that = II 0 . For an acoustically treated 
duct a metric for performance of the lining is the transmission loss defined by 

TL= 201og 10 ^ (93) 

u o 

The presence of possible mode interaction power somewhat complicates the traditional definition 
of transmission loss, and there could also be an argument for replacing II 0 with only the first term 
in equation (90), this being the incident power. 

CONCLUSION 

A reverse flow theorem for acoustic propagation in compressible potential flow has been 
obtained directly from the field equations without recourse to energy conservation arguments. A 
reciprocity theorem for the scattering matrix for propagation of acoustic modes in a duct with either 
hard walls, or a section of locally reacting absorbing wall imbedded in an otherwise hard wall, 
follows. It is found that for a source at a specific end of the duct, suitably scaled reflection matrices 
in direct and reverse flow have a reciprocal relationship. Scaled transmission matrices obtained for 
direct flow and reversed flow with simultaneous switching of source location from one end to the 
other also have a reciprocal relationship. 

The approach presented here is an alternative to the approach of Moehring [3,16], with the 
distinction that no energy conservation condition is used. It has been exploited to provide explicit 
reciprocal relations which are of theoretical interest, but which also have the more pragmatic 
significance of providing a convenient means for benchmarking large scale propagation codes such 
as those used in this investigation. 

Numerical verification of the reciprocal relationships is the subject of a companion paper. 
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Figure 1. A general duct configuration showing volume and surface important in the development. 




RECIPROCITY AND ACOUSTIC POWER IN ONE DIMENSIONAL 
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ABSTRACT 

A reverse flow theorem for one dimensional acoustic propagation in compressible 
potential flow has been obtained directly from the field equations without recourse to energy 
conservation arguments. Reciprocity relationships for the scattering coefficients for propagation 
are derived. It is found that for a source at a specific end of the duct, suitably scaled reflection 
coefficients in direct and reverse flow have a reciprocal relationship. Scaled transmission 
coefficients obtained for direct flow and reversed flow with simultaneous switching of source 
location from one end to the other also have a reciprocal relationship. Reciprocal relations and 
power conservation arguments are used to show that scaled power reflection and transmission 
coefficients are invariant to flow reversal and switching of source location from one end of the 
duct to the other. Numerical verification of the reciprocal relationships is given in a companion 
paper in which multiple mode propagation and one dimensional propagation are considered 

INTRODUCTION 

In a companion paper [1] an acoustic reciprocity theorem in a compressible potential 
flow in non-uniform ducts has been obtained directly from the field equations. It was shown that 
reciprocity holds for ducts with rigid walls and for ducts with a finite length dissipative lining 
imbedded in an otherwise rigid wall. In the present investigation the methodology of [1] is 
applied to the simpler case of one dimensional propagation, valid for low frequencies, that is 
large wave length to geometric cross section dimension ratio. The principal simplification is that 
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acoustic treatment is no longer a possibility within the acoustic formulation. The energy 
formulation of Moehring [2,3] then becomes a straight forward option because dissipative linings 
and cut-off modes do not need to be considered. However, here the approach of [1] based directly 
on the field equations is used . 

Aside from finding a reciprocity relationship, a goal here is to expand upon an invariance 
property for acoustic power transmission in converging-diverging ducts found by Davis [4], This 
turns out to be a result of reciprocity and energy conservation. 

In another companion paper [5], the reciprocity relations found in [1] and in the present 
paper are substantiated by numerical experiments based on a finite element simulation. 


RECIPROCITY IN ONE-DIMENSIONAL FLOW 


In some cases it is possible to use a one-dimensional approximation for steady flow and 
acoustic perturbations. A reverse flow reciprocity relationship can be obtained in this case using 
the method of [1]. Figure 1 shows the duct under consideration, which is in fact three 
dimensional, but it is assumed that the large wave length limit exists so that acoustic propagation 
is one dimensional at each cross section. Furthermore, the duct shape is restricted so that the 
mean flow at a cross section can also be treated as one dimensional. The cross sectional shape of 
the duct is not a consideration in this one dimensional approximation. The one-dimensional 
acoustic continuity equation is 





= 0 


( 1 ) 


where A is the local cross sectional area of the duct. The one-dimensional acoustic momentum 
equation is 


P = 



( 2 ) 
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or 


P = 



(3) 


The steady flow is obtained from 


d_ 

dx 


( P /f) = ° 


(4) 


and the relations between mean flow velocity and mean flow potential and acoustic particle 
velocity and acoustic potential are 


M = 


34 >, 

ax 


(5) 


u 



( 6 ) 


The acoustic relationship between pressure and density is 


P = C r P 


(7) 


Mean flow speed of sound is determined by 


c, 2 = 1 - - A/_ 2 ] 


( 8 ) 


and the mean flow density is 
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( 9 ) 


p, -n -■^ r ^(v4>/v<t> r -M. 2 )j>- 1 

Equations (1) and (9) are non-dimensional using the density and the speed of sound at 
some point, in this case the plane x = 0 . Stagnation conditions could also serve as the 
reference. A reference length R is a (hypothetical) radius corresponding to a (circular) source 
“plane” at x = 0 , the reference cross sectional area. Acoustic potential is non-dimensional with 

respect to c^R, and the acoustic pressure with respect to p jc\ Lengths are made non- 

dimensional with respect to R. Time is scaled with i?/c M - hi the case of harmonic time 
dependence this leads to the definition of non-dimensional frequency r\ r = (oR/c^. co is the 
dimensional source frequency. is the Mach number at the reference point. With this 
convention, the reference non-dimensional frequency and the local non-dimensional 
frequency T| are defined as in [1] where the restriction of one dimensional propagation is not 
invoked. Non-dimensionalization described here differs from the usual one dimensional 
development in which some characteristic duct length is generally used as a reference length. The 
convention used here makes everything consistent with [1], in particular the non-dimensional 
frequency r| . 

Consider the duct shown in Figure 2, which is locally of arbitrary cross section, and is 
terminated on each end by a section of uniform duct, of length sufficient to assure uniform mean 
flow and propagation in terms of acoustic modes. The source plane x = 0 is where the acoustic 
source is specified and the exit plane x - L terminates the duct. For computations the exit 
plane will be assumed non-reflecting. In general, a typical computational problem would seek to 
specify the acoustic field within the duct and the scattering matrix at the source plane for incident 
acoustic modes. Figure 2, originally presented in support of the discussion in [1], can be used 
here to highlight the one dimensional approximation. It is assumed that the mean flow is 
everywhere described in terms of the axial coordinate x so that in the nonuniform section 
M(x,r) = M(x) . Acoustic propagation is assumed to also depend only on the axial coordinate, 
so that it is locally planar. The duct walls do not provide a boundary condition in the one- 
dimensional approximation so that the unit normal n there plays no role. The unit normal on the 
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terminating planes is axially directed. Equations (1) through (9) specify the acoustic field in the 
duct of Figure 2 subject to the restrictions noted. 

Let 4^ e ,r] ‘‘ be a harmonic solution for the acoustic velocity potential for the case of a 
mean flow specified everywhere in the duct by its reference Mach number A? = V<|> r and with 
specified boundary conditions. Let <J> 2 e be a second harmonic solution for exactly the same 
duct with different source conditions, but with the flow reversed, -M = - V^. It is important 
to note that in reversed flow the reference density p r and reference speed of sound c r are 
unaltered. Because of equation (1), in the case of a harmonic source at non-dimensional 
frequency rj , it follows that 

f {4> 2 [*VPi + V-CpAV^ + MA p t )] - <b l [ir\Ap 2 * V-(p/V(j> 2 - MAp 2 )]}dx = 0 

( 10 ) 

Here, for ease of notation, and for a degree of similarity with [1], the gradient operator is used to 
denote 



and 


( 11 ) 


M = M e (12) 

r r x 

with e the unit vector in the x direction. With application integration by parts, which is 

X 

equivalent to the divergence theorem in [1], and by using equation (2) to partially replace p ( 
and p 2 , and to subsequently eliminate the remaining integral on x , equation (12) is reduced 
to 


3<|> 


3cL 


[ ^ iP ' A -d + ^ Pl) ' 3 ° = 0 


(13) 
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When this is written entirely in terms of acoustic potential the result is 


[Ap {<J>,[(1 - M 2 )-0 - - 4>! [( 1 - AT 2 )— + ztiA/4> 2 ]} ] q - 0 (14) 

In [1], eigenfunction expansions are used to define the acoustic potential field in the uniform 
flow regions at x — 0 and x - L . This is still appropriate, with the simplification that only the 
plane waves propagate, and no analysis is required to determine the modes and wave numbers. 

There is one “mode” propagating in each direction in the uniform duct sections, so that 
there is only a superposition of a right and left wave, given for example here for the nominal 
flow solution at x = 0 in the form [1] 


(b = a/ + a, 


(15) 



(-ik~ ) a 

V X i 1 


(16) 


The axial wave numbers are 


/ 



n ) 


— - — -( - M± 1) 
1 - M 2 


(17) 


/ 



V 


— r( 1) 

1 - M 2 


(18) 


The superscript choice ± corresponds to right and left waves. r\ and M represent local values 
of non-dimensional frequency and Mach number determined by the local speed of sound, with 
q still based on the reference length. The operator arising in equation (14) 
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!{<{)} - (1 - M 2 )^- - ir\M$ = a<]> 


(19) 


introduces the parameters a* and a' corresponding to right and left waves, for example at 
x = 0, 


a \ = P r [-/(1 ~ ~ ZT| ^ = ( 20 ) 

v-\ = P r [-/(1 ~ M 2 )k~ - ix\M ] = /p r ri (21) 

Similar expressions exist for a~ , in reversed flow. It is noted that a 2 = a i =a 

a" = a" = - a . a is evaluated at the ends x = 0 and x = L as required. 

2 1 

With these definitions, equation (14) can be rearranged in a matrix form [1] , 


«r 

T 

A o% 


a 2 


a 2 _ 

T 

o 

o 




> 


< 


= < 


> 



> 



A L a L 




K , 


A i a L 


A. 


( 22 ) 


Modal amplitudes a* and a~ and b* and b~ , acoustic potential wave amplitudes at 
x = 0 and x = L , are related by the acoustic potential scattering matrix according to 




(23) 


where the scattering matrix is defined as 


[ 5 ] 


R T 
T R 


(24) 


7 


158 


R and T are the reflection and transmission coefficients for a source at x - 0 with a 
reflection free termination at x - L . R and T are reflection and transmission coefficients 
for a source at x = L and reflection free at x = 0 . There will be a scattering matrix [S^ ] for 
nominal mean flow and a second one [&,] for reversed flow. Equations (23) and (24) can be 
used in equation (22) to obtain 



T 

4- 


+• 

T 


< 

■ [S,] r M[cc]' 

a 2 

> = < 

a 2 

> [S 2 f[^][a]< 

a i 

► 


b 2 \ 


S , 


b l\ 


where the definition 




A L a L 


[A][a] 


(26) 


is introduced. The scattering matrices for nominal and reversed flow are [S) ] and [>S 2 ] . The 
implications of equation (26) are summarized by the reciprocity relations 

([5 l fM][a]) = [A][a][S 2 ] (27) 


and 


[A] [a] [SJ = ( [^][«][5 2 ]) r 


(28) 


Equation (27) and (28) are used to establish the following relationships for the acoustic 
potential reflection and transmission coefficients: 


= * 2 (29) 

R { = R 2 ( 3 °) 
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(31) 

(32) 


^ 2^0 a 0 ^l^L a L 

f A% ~ ^2 ^L a L 


By using the definitions of a Q and a £ Defined by equations (20) and (21) evaluated at x-0 
and x = L , equations (31) and (32) are written 


A pc 
T _ _0 _£ 

1 ^4, p c 

L r r n 


v Pr v 


j p c 
T - _£ 

2 d, P C 

L % r o 


' P ) 
_2 

V % 


(33) 


(34) 


Equations (29 ) and (30) produce the interesting result that the acoustic potential reflection . 
coefficients, at the left end and at the right end are invariant to flow direction. No reciprocal 
relationships link the left to right transmission coefficients 7^ and T 2 or the right to left 
transmission coefficients T ' and in nominal flow and reversed flow. Equation (33) 

link s the left to right transmission coefficient T { in nominal flow to the right to left 
transmission coefficient T 2 in reversed flow. Equation (34) links the left to right transmission 
coefficient in reversed flow T to the right to left transmission coefficient T ^ in nominal flow. 

The reciprocal results of equations (29 )-(32) are for acoustic potential modal amplitudes. 
In the case of acoustic pressure modal amplitudes a transition to acoustic potential modal 
amplitudes is made through equation (3). The transformations from acoustic potential to acoustic 
pressure for right and left waves in nomnal and reversed flow are 

<b\ = P *Pi > = P ~pi 

<J )’ 2 = P ’ ^"2 = P ^ Pi (36 ^ 
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where 


P* = 


1 + M 
~ ir )P r c r 


P' = 


1 - M 

-^7r 


(37) 


Mach number and non-dimensional frequency are based on the local speed of sound. Non- 
dimensional frequency is still based on the reference length. Non-dimensional density and non- 
dimensional speed of sound are evaluated locally. 

Equations (35) and (36) are used to introduce pressure modal amplitudes in equation (22). 
The scattering matrix is now in terms of pressure scattering coefficients 


[5] 


R T 
T 1 


(38) 


After following the same steps leading to equations (29)-(32), reciprocity relations for acoustic 
pressure scattering coefficients are found to be 


p- _ p* 
Y. A a — - = R A a — - 

100n+ 2UUr>- 

P 0 ^ o 


or 


- a + M J ^ 

R , = — R , 


1 ( 1 _ M ) 


2 2 


(39) 


_ py - p'l 

R, A. a = R A a 

1 L L r, P‘, 


or 


- (i-^r T 

R = — R 

1 (1 + M r ) 2 2 


(40) 


py - p 0 

T.A. a r — = 77 A. a 

1 “ f, 2 0 " p-. 


' a l P, c ,(l*W t ) 2 2 

n n 


(41) 
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p- _ P'i 

^\%^- T , A L a L — 


or 


P‘ 


P' 


_ A P r c r (1 - M ) 

TjT _ U o o U 'J' 

2 An/*/ i 1/ \2 1 


h P r c r (1 - MY 


(42) 


Equations (39)-(42) are written both implicitly in terms of coefficients defined above and 
explicitly in terms of local area and flow conditions. The mach number is for the nominal flow. 

Pressure reflection coefficients are not invariant in reversed flow but are reciprocally 
related as given by equations (39 ) and (40). No reciprocal relationships link the left to right 
transmission coefficients and Y 2 or the right to left transmission coefficients f, and 

Y 2 in nominal flow and reversed flow. Equation (41) links the left to right transmission 
coefficient Y in nominal flow to the right to left transmission coefficient ? 2 in reversed flow. 
Equation (42) links the left to right transmission coefficient in reversed flow T 2 to the right to 
left transmission coefficient T } in nominal flow. 

ACOUSTIC POWER CONSIDERATIONS 

In the case of one dimensional propagation, additional results can be obtained by the 
consideration of acoustic power. The definition of acoustic intensity due to Morfey [6] is valid as 
part of a conservation law in non-uniform ducts for compressible potential flow. This definition 
of acoustic intensity is used in the uniform flow sections on either end of the nonumformity to 
obtain acoustic power expressions. 

For propagation in uniform flow the Morfey intensity formulation simplifies to the scalar 

form 



( 1 + M 1 ) pu + p r c r Mu 2 + 


P c 

r r r 


Mp- 


Acoustic power is obtained by integration over a cross section, to yield 
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— = f f I (l + M 2 ) pu + p c Mu 2 + Mp 2 \ dS (45) 

A re/P- C ~ Vi ^ I 

A modal expansion as given by equations (26) and (27) is used to obtain an acoustic power in 
terms of acoustic potential modal amplitudes. The result is 


= <T />,?«* +a-'P~ u a- (46) 

3 U 11 

A ,p c 

where, for example, cl denotes the complex conjugate of u . Power can be represented in 
terms of acoustic potential amplitudes as in equation (46) in uniform duct sections at x - 0 and 
x -L, The power coefficients P ** and are one dimensional analogs of power matrices 

introduced in [1]. No acoustic power is attributed to the interaction of acoustic modes because 
there are no cut-off modes in this one dimensional model. The power coefficients can be easily 
obtained by using the results for the plane wave mode from the general expressions of [1], 


it 


p, c , 


11 


— OC TL 2 A. 

o r r r '0 0 


(47) 


Energy conservation arguments lead to the conclusion that II. ^ This is 

the traditional result, simplified from the result of [1] because of the absence of interaction power 
associated with cut-off modes. Incident power at x = 0 is given by 


n 

inc 


= fl *' r:p;: R,a* 


a 


Reflected power in nominal and reversed flow at x ~ 0 is given by 


Kf- a " R : p u/y 


r ; p n 0 V 


(48) 
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Transmitted power in nominal and reversed flow x = L is given by 


n 


trans 


* t: p:: r,a* 


i a 


n = a * t: p:: 

trans l 1 1 




(48) 


Subscripts on the power coefficients denote the uniform section in which they are evaluated. 
From equation (29) R x = i? 2 , that is, the acoustic potential reflection coefficient is invariant 
to flow direction, as is P . Therefore, the power reflection coefficient defined by either 


R 


R \ P n Pj 

n 


ip. 


or 


R 


R ; P n R 2 

n 


p ' l ' "B p„ 

11. U n 


I*, 


.(49) 


The power reflection coefficient is invariant to flow direction, P^ = P^. It is concluded that the 
power transmission coefficient are also invariant to flow direction, T = T^. A similar 
development based on equation (30) would show that P_ = P^ and therefore that ^n 0 ' 

It is now possible to connect the power transmission coefficients when the source is 
moved from one end of the duct to the other. For a source at x = 0 the power transmission 
coefficient (the ratio of transmitted to incident power) is in nominal flow 




Pr 

1 


A 

c f |r ' 

C r A 0 


(50) 


For a source at x - L in reverse flow the power transmission coefficient is defined for right to 
left power according to 


T 

TC 


0 


T m P T 
2 11 2 




u 

13 


p c A 

0 | f | 2 
p c A n 

~ r r 0 


( 51 ) 
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Equations (50) and (51) are deduced by using the definitions of incident and transmitted power 
appropriately evaluated at x = 0 or x = L. The reciprocal relationship of equation (3 1) is used 
to replace f in equation (51) yielding the result 


P r ^ A L 2 

r = — — -r l ^ I 

R 0 P C A. 1 

r n r u 0 


(52) 


= f = f and 


It has thus been shown that T = T . It is therefore concluded that T - T 

*9 *9 _ „ 9 9 

thence from power conservation arguments that R = R = R - R ■ This completes the 

r K e Tig n 9 Tig 

interesting result that for the one dimensional model power reflection and transmission 


coefficients are invariant to flow reversal and switching of the source location. 

The result that T = f states that the transmission coefficient for the nominal flow 

*9 *9 

direction is the same from either end of the duct. This result contains the invariance theorem of 
Davis [4], but also admits a generalization. Davis found that for a converging-diverging non- 
uniformity in an otherwise uniform duct, for equal pressure amplitude input at the upstream and 
downstream ends of the duct, the transmitted power is related by 


_ n +M ) 2 (53) 

flj (1 -Mf 

Here II 1 is the transmitted power at the downstream end due to a source at the upstream end 
and fl 1 is the transmitted power at the upstream end due to a source at the downstream end, 
both in the nominal flow. Because the result of Davis is for a converging diverging duct, with 

both ends of the same area, M Q = M L = M. 

By using equations (46) and (47), which define acoustic power in terms of acoustic 
potential amplitudes, and by modifying them using equations (35)-(37) to relate acoustic pressure 
amplitude to acoustic potential amplitude, the incident power at x = 0 is 

n ,J* = o) = Ib;i 2 -^(i*M 0 )2 ( 54 ) 

r r r 

a n 

Similarly, for an acoustic pressure input at x = L , the incident power there is 
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(55) 


n. (x = l) = hpr '\ 2 — - 

mc v 2 L 0 C 

r r r 


(l 



p Q * and p L are pressure mode amplitudes incident at the ends x - 0 and x - L. 
Superscripts + and - reinforce the idea that at x = 0 the incident mode is a right running wave 
and at x - L the incident mode is a left running wave. The tilde ( - ) reminds that the source is 
at the end x = L. By using the fact that the transmitted power in each case is the product of the 
incident power and the appropriate power transmission coefficient, which is the same for either 
case ( T = T ) and that the incident modal pressure amplitudes are the same for either 

v *e *6 

source ( [p 0 *| = \p L | ) > it follows that 


n 

trans 


n 

trans 


-4, \ % uy 

A L P r C r (1 -M L f 
ti n L 


.(56) 


' Equation (56) contains Davis’ result [4] in the case of a converging diverging duct when the duct 
area, flow density, speed of sound and Mach number are the same at both ends. Other results are 
possible involving acoustic potential and acoustic pressure amplitudes from the core result that 
the power transmission coefficients are invariant. 

CONCLUSION 

A reverse flow theorem for acoustic propagation in one dimensional compressible 
potential flow has been obtained directly from the field equations without recourse to energy 
conservation arguments. A reciprocity theorem for the scattering coefficients for propagation of 
acoustic modes has been obtained. Reciprocal relations and power conservation arguments are 
used to show that scaled power reflection and transmission coefficients are invariant to flow 
reversal and switching of source location from one end of the duct to the other. 

Numerical verification of the reciprocal relationships has been made using a finite 
element model for duct propagation, and is reported in a companion paper [5]. 
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one dimensional approximation, as is the outward unit normal on the duct walls. 
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ABSTRACT 

A reciprocity theorem for the scattering matrix for propagation of acoustic modes in a 
duct with acoustically hard walls or with acoustically absorbing walls has been given in a 
companion publication. It was found that for a source at a specified end of the duct, suitably 
scaled reflection matrices in direct and reverse flow have a reciprocal relationship. Scaled 
transmission matrices obtained for direct flow and reversed flow with simultaneous switching of 
source location from one end to the other also have a reciprocal relationship. A reverse flow 
theorem for the equivalent one dimensional propagation model, which is a good approximation 
to the three dim ensional model at low frequencies, was also obtained. In this case, using 
reciprocity and acoustic power conservation arguments it is additionally found that the acoustic 
power transmission coefficient is the same for a source at either end of the duct for a given flow 
direction. This result leads to an invariance theorem which relates acoustic power propagated due 
to sources of equal pressure amplitude at the two ends of the duct. Numerical verification of 
these reciprocal relationships is given here for propagation in axially symmetric (circular and 
annular) ducts with multi-modal propagation and at low frequencies when a one dimensional 
model is appropriate. 

INTRODUCTION 

In two companion papers [1,2], a reverse flow reciprocity theorem has been developed for 
acoustic propagation in non-uniform ducts carrying compressible mean flow. In [1] the general 
case of multi-mode propagation is considered and reciprocity is shown for ducts with either rigid 
walls or with walls which include a normally reacting, dissipative section. For a source at one 
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end of the duct, scaled reflection matrices in direct and reverse flow have a reciprocal 
relationship. Scaled transmission matrices obtained for direct flow and reversed flow with 
simultaneous switching of source location from one end to the other also have a reciprocal 
relationship. In [2] reverse flow reciprocity is shown for the long wave length approximation 
when a one dimensional model is applicable. In this case acoustic treatment is not part of the 
model. Results similar to the multi-modal case are established for reciprocal relationships for 
reflection and transmission coefficients. Additional results which are part of a general power 
transmission invariance principal are also found as a result of reciprocity and energy 
conservation. This invariance principal contains as a special case a result found by Davis [3]. 

The reverse flow reciprocity theorem is developed directly from an integral relationship 
based on the acoustic field equations, using an approach similar to that used in [4] in the case of 
propagation in a non-uniform duct in the absence of flow. In particular, the development does not 
begin with an energy principle. This is in distinction to the approach of Moehring [5,6]. The 
major complication which arises in the present formulation is the case when a portion of the duct 
wall is acoustically treated with a normally reacting dissipative lining. The boundary condition of 
Myers [7] is manipulated, using identities of vector calculus suggested by Moehring [6], to make 
it possible to establish reciprocity in this case. 

The reverse flow reciprocity theorem is developed in detail in references [1] and [2] and 
the results are briefly summarized here. Figure 1 shows a non-uniform duct section bounded on 
either end by uniform sections (long enough to have essentially uniform flow so that acoustic 
propagation can be synthesized in terms of duct modes). At the two ends of the duct the acoustic 
field is the superposition of modes propagating to the right and to the left (including cut off 
modes which technically do not propagate, but which can be segregated into right and left 
modes). Amplitudes a* and a~ refer to right and left modes at the end x = 0 and b n 
and b ‘ refer to right and left modes at x = L . a* , a' ,b* ,b ,are vectors of modal 

n 

amplitudes. These modal amplitudes are related by the scattering matrix [S'] according to 
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The scattering matrix is defined as 


[S] 


[R] [f] 
[T] [R] 


( 2 ) 


Contained in [5] are the usual reflection matrix [R] and transmission matrix [T] for acoustic 
modes incident at x = 0 and reflection and transmission matrices [R] and [f] for modes 
incident at x = L. In multi-modal propagation the scattering matrix relates all modes which are 
considered. In the case of one dimensional propagation (the long wave length approximation), 
the scattering matrix relates only two modal amplitudes at each end. The reflection and 
transmission matrices are scalars, defined as reflection and transmission coefficients. 

In the context of reversed flow reciprocity there is a scattering matrix [S) ] for nominal 
mean flow and a second one [&,] for reversed flow. It is the relationship between [S y ] and 
[5 2 ] which is considered in [1,2]. Modal amplitudes in the present discussion are in terms of 
acoustic potential duct modes, because the acoustic field equations are naturally in terms of 
acoustic potential. Equivalent results are obtained in [1,2] for acoustic pressure modal amplitudes 
and it is only necessary at this point to refer to the properties of the acoustic potential scattering 
matrices. 

The reverse flow reciprocity principle is 


[$,J r [•/][«] = M[a][S 2 ] 


( 3 ) 


or 


[•/][«][$,] - (M[«][S 2 ]) r (4) 

The diagonal matrices [ J] and [a] are scaling matrices which have elements depending on 
the mode considered. The evaluation of these matrices is covered in detail in [1,2]. 

Equations (3) and (4) show that a weighted version of the nominal flow acoustic potential 
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scattering matrix and similarly weighted version of the reversed flow acoustic potential scattering 
matrix are transposes of one another. In terms of the acoustic potential reflection and 
transmission coefficient matrices the result is 


K] E* 2 ] 

(5) 


(6) 

[r l ] r [J t ][« t ] 

(7) 

[f,] r [J- 0 ][«„] -V L ] K][r 2 ] 

(8) 


Subscripts 0 and L refer to the evaluation of the relevant scaling coefficients at the two ends 
x = 0 and x = L . The reciprocal relationships of equations (5)-(8) involve acoustic potential 
reflection and transmission coefficient matrices, with diagonal elements representing reflection 
and transmission coefficients in the incident modes (here referred to as direct reflection or 
transmission) and off diagonal reflection and transmission coefficients from the incident mode to 
another mode. Equations (5) and (6) show that direct acoustic potential reflection coefficients 
are invariant in reversed flow. The transmission coefficient matrix pairs [T x ] , [ T 2 ] and 
[f] , [? 2 ] are not reciprocally related but the pairs [TJ , [f 2 ] and [f ( ] , IT 2 ] are related 
by equations (7) and (8). The notation convention uses the subscripts 1 and 2 to denote flow 
direction (1 being nominal, 2 being reversed). Tilde, or lack thereof denotes source location (tilde 
denoting source location reversal). 

In the one dimensional approximation [2] the same reciprocal relations apply, but 
in a simplified form. Reflection coefficients in nominal and reversed flow are invariant: 
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( 9 ) 


R { = R 2 (10) 

Scaled transmission coefficients are invariant to simultaneous flow reversal and source plane 
reversal: 




( 11 ) 


P r 


Pr. 


TA r — = T,A n — 


2 L 


10 


( 12 ) 


In the case of one dimensional propagation there are additional results which can be 
deduced based on reciprocity and energy conservation. Power transmission coefficients T n 
are defined as the ratio of transmitted power to incident power. It is found that 


T = T = f = f (12) 

", "2 "2 

Here, as in the previous discussion the subscripts 1 and 2 refer to flow direction and the tilde or 
lack thereof refers to source location. Power transmission coefficients are invariant to flow 
reversal and source location reversal. That is, the power transmission coefficient is the same for 
flow in either direction, for a source at either end of the duct. Power reflection coefficients R ^ 
are defined as the ratio of reflected power to incident power. It is also found that 


R 


7T 


1 




(14) 


Power reflection coefficients are invariant to flow reversal and source location reversal. 
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The results summarized by equations (5)-(14) are interesting theoretically and also 
provide useful benchmarks which can be used to validate propagation calculations. In the 
following sections two finite element codes for duct propagation, one multi-modal and the other 
one dimensional, are used to demonstrate several of these reciprocal relations. 

ACOUSTIC PROPAGATION IN A COMPRESSIBLE POTENTIAL FLOW 

Reciprocity relations previously described will be verified by computations based on two 
FEM codes for duct propagation, one for multi-modal propagation and the other specialized for 
one-dimensional propagation. In this section only a brief description of the multi-modal 
propagation code will be given. Details of the FEM modeling approach can be found in 
references [8,9]. 

A formulation in terms of acoustic potential is used to produce a weak formulation 
suitable for finite element discretization to reduce the field equations to a single scalar variable. 
The geometry of the duct in Figure 1, and the steady flow field is axially symmetric. The 
acoustic field is not axially symmetric but is represented as azimuthally periodic in a cylindrical 
coordinate system with x being the axis of symmetry, r the cylindrical radius in a circular 
cross section at x = 0, and 0 the angular coordinate. Solutions are sought in angular harmonics 
of a Fourier Series in 0 enumerated by the angular mode number m . This reduces the solution 
domain to a two dimensional x , r plane, shown in Figure 1. The duct shape in a 0 = constant 
plane shows the surface S which defines the duct shape and could include an inner surface for 
an annular duct. Part of S includes S , which is a locally reacting acoustic treatment. 

The acoustic field is assumed to be harmonic in time at non-dimensional frequency r^. 
Geometry is non-dimensional based on a reference length generally chosen as the radius of the 
inlet at the source plane, R. Acoustic and steady flow variables are non-dimensional based on 
reference values of the speed of sound and density of the medium, , c_, generally defined at 
the plane of the acoustic source. The non-dimensional frequency is r\ r = oo R / , with (othe 

harmonic source frequency. 

Field equations for continuity and momentum and the isentropic equation of state are 
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used in a weighted residual statement to obtain an integral formulation which is then written in 
discrete form using standard FEM procedures. In terms of acoustic potential the weak 
formulation is 



+ ir\ r [W(V r -V <j>) -(F-VFF) <J>] - r\ 2 W$}dV 


P , 


= jf-L{c 2 WV$ - V r W(V r -V< J)) - ii\ r VW$}-ndS 


(15) 


where the local non-dimensional steady flow velocity is F = V (J>^ , with (j) r the non- 
dimensional steady flow velocity potential. The local non-dimensional density and speed of . 
sound are p , c . The surface integral on the right hand side introduces the noise source and 
termination conditions on S Q or S L and a possible impedance boundary condition on S inside 
the duct. In the present investigation it is the impedance boundary condition which is of interest 
on S w , a. portion of S. In equation (15), the weighted residuals statement, W represents an 
arbitrary weighting function selected from the class of continuous functions. In this weak 
formulation the approximation to the solution (J) is also chosen from the class of continuous 
functions 

At a duct wall the mean flow is tangential to the wall and F_ • « = 0 causing the 
boundary integral (the contribution to the right hand side of equation (1) related to the impedance 
condition) to become 


I b = f f p r fVV<p -ndS 


(16) 


It is shown in [10 ] that at a wall of admittance A the weighted residual boundary integral of 
equation (16) on the duct surface^ , derived from the Myers boundary condition [7], is 
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An accurate representation of this impedance boundary condition is essential to obtaining 
verification of reciprocal relationships when acoustic treatment is inserted in the duct wall [10]. 
References [8,9 ] deal with propagation and radiation to the far field from open ended ducts. In 
the code described here, rather than model radiation to the far field from the open end, imposes 
non-reflecting boundary condition at the termination. The boundary integral of equation (15) is 
used to introduce the source (at either end of the duct) as a superposition of acoustic potential 
duct modes and to implement the non-reflecting boundary condition based on another 
superposition of duct modes. The one-dimensional code is based on the one-dimensional field 
equations [2], and therefore has no provision for acoustic treatment. Other details of this code are 
similar to the multi-modal code. In both cases the steady mean flow field which is required data 
for propagation calculations is provided by an FEM potential flow code which introduces 
compressibility by iteration of successive incompressible flow problems. Steady flow is 
produced on the acoustic FEM mesh for convenience of data transfer. 

The FEM duct codes provide solutions for the acoustic potential field which is post- 
processed to obtain acoustic pressure. Included as part of the solution are the acoustic potential 
modal amplitudes a' , a~ ,b* ,b~ . These are also converted to acoustic pressure modal 
amplitudes. Additional post-processing produces computations of acoustic power and power 
reflection and transmission coefficients. 

NUMERICAL EXPERIMENTS ON RECIPROCITY 

The first numerical verification of the reciprocal characteristics of the scattering matrix 
for an axially symmetric duct has been carried out in the case of a duct with a transition from 
annular to circular, as shown in an x , r slice in Figure 2. The interior contour is that of a typical 
turbofan inlet and the uniform extensions are added to meet the requirements of the present 
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analysis (uniform flow and proper definition of acoustic eigenfunctions). The finite element 
mesh used in the computations is shown on this figure and is typical for examples cited here. The 
conditions for the “nominal” case are standard atmospheric conditions at the source plane (x = 0 
), Mach number at the source plane = 0.27 , directed left to right (opposite to the direction 
in an inlet). The non-dimensional frequency (based on the source plane) is T| r = 10 . In the 
nominal case the input plane for scattering is the source plane at x = 0 . Figure 3 shows iso- 
potential contours for the steady compressible flow in the duct which varies roughly between 
M = 0.27 and M = 0.15. Figure 3 is unaltered in form for reverse flow (right to left and in the 
direction expected in an inlet). The acoustic analysis is based on input modes with angular 
dependence m = 3 , for which there exist two propagating modes at each end of the duct (n=l, 
n=2).The third mode, n=3, is cut off at both ends of the duct with cut-off ratio £ = 0.87 at 
x = 0 and £ = 0.84 at x - L . Two cases of duct wall characteristics are studied. In the 
first case the duct walls are acoustically hard, that is the impedance is infinite and the admittance 
vanishes. In the second case the outer duct wall is acoustically treated from x = 1.0 through 
80 % of the non-uniform section. The impedance is chosen as Z = 2.0 - 1.0 i , which is not 
optimum for attenuation for the given conditions, but is not untypical for aircraft applications. 
The acoustic power attenuation with the simplest radial mode incident (m— 3, n=l) is about 9 dB, 
so there is a significant decrease in acoustic power from one end of the duct to the other, 
attributable to the wall treatment. 

To generate the scaled reflection and transmission coefficient matrices in this case input 
modes n - 1 , n — 2 and n - 3 are considered. This produces 3x3 matrices which include 
two propagating modes and one mode which is cut off at both ends. To build the reflection and 
transmission matrices to verify equations (5) and (7) it is required to consider nine propagation 
cases: three input modes at x - 0 for nominal flow; three input modes at x = 0 for reverse 
flow; and three input modes at x = L for reverse flow. 

Reference to contours of equal acoustic pressure provide evidence of the differences 
induced by varying the source mode number, source location and the flow direction. Figure 4 
shows acoustic iso-pressure contours for the duct with acoustic treatment when the input radial 
mode at x - 0 is n = 1 . There is significant scattering and this is verified by reference to the 
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scattering coefficients. Figure 5 shows the case of mode n = 3 input at x - L in reverse 
flow. This mode is cut off at both ends, and it should be noted how rapidly the acoustic pressure 
level is attenuated away from the source plane. 

Table 1 shows the scaled acoustic potential reflection coefficient matrix in nominal flow 
with the scattering plane (acoustic source plane) at x = 0 , which corresponds to the left hand 
side of equation (5). Table 2 is the reflection matrix at x = 0 in reversed flow, corresponding to 
the right hand side of equation (5). Equation (5) predicts that these matrices should be 
reciprocals of one another. Tables 1 and 2 verify this with exceptional accuracy. 


mode 

1 

2 

3 

1 

-0.14380 - 0. 03098 i 

-0.13986 + 0.07217 i 

-0.153(-4) - 0.549(-5) i 

2 

-0.13746 + 0.71323 i 

-0.04113 + 0.20572 i 

0.104(-3) + 0.732(-5)i 

3 

-0.424(-5) + 0.841(-5) i 

-0.144(-4) - 0.903(-4) i 

-0.682(-5) + 0.497(-8) i 


Table 1 . Scaled acoustic potential reflection coefficients for nominal flow and source at x = 0 
for the transition from an annular to circular hard wall duct. Corresponds to the left hand side of 
equation (5). r\ = 10 , m = 3. 


mode 

1 

2 

3 

1 

-0.14380 - 0. 03098 i 

-0.13746 + 0.71323 i 

-0.424(-5) + 0.841(-5) i 

2 

-0.13986 + 0.07217 i 

-0.04113 + 0.20572 i 

-0.144(-4) - 0.903(-8) i 

3 

-0.153(-4) - 0.549(-5) i 

0.104(-3)+ 0.732(-5) i 

-0.682(-5) + 0.497(-8) i 


Table 2. Scaled acoustic potential reflection coefficients for reverse flow and source at x = 0 
for the transition from an annular to circular hard wall duct. Corresponds to the right hand side of 
equation (5). v\ r = 10 , m = 3. 

Tables 3 and 4 verify the prediction of equation (7). Table 3 gives the scaled acoustic 
potential transmission coefficients in nominal flow with the source at x - 0 . Table 4 gives the 
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scaled transmission coefficients in reversed flow with the source shifted to x - L . Equation (7) 
predicts a reciprocal relationship which is accurately substantiated in Tables 3 and 4. A point of 
interest in the results shown in Tables 1 through 4 is that reciprocity extends to cutoff modes in 
which case the power transmission is only accounted for by interaction of left and right modes. 


mode 

1 

2 

3 

1 

-4.12246 + 0.00532 i 

-0.56880 + 0.23281 i 

0. 1 5 1 (-4) - 0.109(-4) i 

2 

0.58902 + 0.37376 i 

-1.74565 - 0.27580 i 

0.834(-4) + 0. 1 74(-4) i 

3 

-0.300(-4) - 0.897(-4) i 

-0.169(-3)+ 0.105C-3) i 

0.808(-8) - 0.488(-8) i 


Table 3. Scaled acoustic potential transmission coefficients for nominal flow and source at 
x - 0 for the transition from an annular to circular hard wall duct. Corresponds to the left hand 
side of equation (7). = 10 , m - 3. 


mode 

1 

2 

3 

1 

-4.12246 + 0.00532 i 

0.58902 + 0.37376 i 

-0.300(-4) - 0.897(-4) i 

2 

-0.56880 + 0.23821 i 

-1.74565 - 0.27580 i 

-0.169(-3) - 0.105(-3) i 

3 

0.15 1(-4) - 0.109(-4) i 

0.834(-4) + 0.174(-4) i 

0.808(-8) - 0.488(-8) i 


Table 4. Scaled acoustic potential transmission coefficients for reverse flow and source at x = L 
for the transition from an annular to circular duct. Corresponds to the right hand side of equation 
(7). t| r = 10 , m - 3 . 

Tables 5 through 8 are for the case with acoustic treatment in place. Table 5 shows the 
scaled acoustic potential reflection coefficient matrix in nominal flow with the scattering plane 
(acoustic source plane) at x = 0 , and with the nonuniform portion of the duct outer wall 
acoustically treated. This corresponds to the left hand side of equation (5). Table 6 is the 
reflection matrix at x = 0 in reversed flow, corresponding to the right hand side of equation (5) 
for the same acoustic treatment configuration. Equation (5) predicts that these matrices should 

11 


180 





be reciprocals of one another. Tables 5 and 6 verify this, again with exceptional accuracy 


mode 

1 

2 

3 

1 

0.16904 - 0.23900 i 

-0.12975 - 0.18014 i 

-0.482(-3) - 0.677(-4) i 

2 

-0.07368 - 0.10133 i 

-0.00984 - 0.03495 i 

-0.275(-3) + 0.950(-4) i 

3 

-0.305(-3) + 0.632(-5) i 

-0.283(-3) + 0.556(-5) i 

-0.693(-5) + 0.51 0(-6) i 


Table 5. Scaled acoustic potential reflection coefficients for direct flow and source at x = 0 for 
the transition from an annular to circular acoustically treated duct. Corresponds to the left hand 
side of equation (5). t| = 10 , m - 3 . 


mode 

1 

2 

3 

1 

0.16904 - 0.23900 i 

-0.07368 - 0.10133 i 

-0.305(-3) + 0.632(-5) i 

2 

-0.12975 - 0.18014 i 

-0.00984 - 0.03495 i 

-0.283(-3) + 0.556(-5) i 

3 

-0.482(-3) - 0.677(-4) i 

-0.275(-3) + 0.950(-4) i 

-0.693(-5) + 0.510(-6)i 


Table 6. Scaled acoustic potential reflection coefficients for reverse flow and source at x = 0 
for the transition from an annular to circular acoustically treated duct. Corresponds to the right 
hand side of equation (5). r| = 10 , m = 3 . 

Tables 7 and 8 verify the prediction of equation (7) in the case of an acoustically treated 
outer wall. Table 7 gives the scaled acoustic potential transmission coefficients in direct flow 
with the source at x - 0 . Table 8 gives the scaled transmission coefficients in reversed flow 
with the source shifted to x = L . Equation (7) predicts a reciprocal relationship which is 
accurately substantiated in Tables 7 and 8. Again, the applicability of the reciprocity principle to 
cut-off modes is verified. 
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mode 

1 

2 

3 

1 

-1.21549 +0.02495 i 

-0.25788 + 0.60030 i 

0.780(-4) - 0.227(-3) i 

2 

0.38577 - 0.37557 i 

-0.76395 + 0.24663 i 

0.146(-3) - 0 ,741(-4) i 

3 

-0.319(-4) - 0.753(-4) i 

-0.532(-4) + 0.849(-4) i 

0.109(-7) - 0.205(-7) i 


Table 7. Scaled acoustic potential transmission coefficients for direct flow and source at x = 0 
for the transition from an annular to circular acoustically treated duct. Corresponds to the left 
hand side of equation (7). r\ r = 10 , m = 3 . 


mode 

1 

2 

3 

1 

-1.21549 + 0.02495 i 

0.38577 - 0.37557 i 

-0.319(-4) - 0.753(-4) i 

2 

-0.25788 + 0.60030 i 

-0.76395 + 0.24663 i 

-0.532(-4) + 0.849(-4) i 

3 

0.780(-4) - 0227(-3) i 

0.146(-3) - 0.741(-4) 

0.109(-7) - 0.205(-7) i- 


Table 8. Scaled acoustic potential transmission coefficients for reverse flow and source at x - L 
for the transition from an annular to circular acoustically treated duct. Corresponds to the left 
hand side of equation (7). = 10 , m = 3 . 


The next case considered involves a steady flow in which the Mach number becomes 
relatively high, emphasizing the dependence of the acoustic treatment boundary condition on 
Mach number. An additional complication introduced here is the segmenting of the acoustic 
treatment into two equal length parts with different impedances, spanning the entire transition 
section. A continuous transition occurs within one element of the FEM mesh. The impedances in 
this case are Z^ = 2.0 - 1.0 i , Z^ = 3.0 - 2.0 i (numbered left to right) . This satisfies 
the requirement of the reciprocity theorem and also simulates a near discontinuity of impedance. 
Here a converging duct with a contraction ratio of o = 0.5 , as shown in Figure 6, accelerates the 
flow (iso-potential contours are shown on Figure 6) from a Mach number at the nominal source 
plane of M = 0.13 to M = 0.71 at the exit plane. An acoustic propagation analysis has been 

13 


182 





carried out for the non-dimensional frequency = 10 for a source with angular mode m = 3. 
In this geometry and the resulting steady flow there are two propagating modes at x = 0 
(determined in the hard wall case), but just one propagating mode at x = L . Results given here 
are for scaled potential reflection and transmission matrices. 3x3 reflection and transmission 
matrices are investigated by considering input radial modes n = 1 , n = 2 , and n = 3 . At 
x = 0 the two propagating modes and one cut-off mode have cut-off ratios 
£ = 2.40 , l = 1.26 , 5 = 0.89, respectively. At x = L the single propagating mode and two 
cut off modes have cut-off ratios § = 1.65 , § = 0.87 , £ = 0.61. An interesting feature of 
this geometry and flow is that mode n = 2 makes a transition from cut on to cut off in going 
from left to right. Mode n = 3 is deeply cut off at x - L . 

Figure 7 shows acoustic iso- pressure contours for the case of nominal flow (left to right) 
with the mode n = 1 input at x = 0 . The contours are consistent with a well cut on mode 
and significant scattering. Figure 8 shows contours for the case of reverse flow (right to left), 
with the source at the end x = L , and the input mode n = 3 . This mode is deeply cut off and 
it should be noted how rapidly the acoustic pressure levels decay away from the source plane. It 
can be concluded that this mode effectively produces no acoustic pressure at x = 0 . 

Tables 9 through 12 are presented to verify the predicted reciprocal characteristics of the 
scaled acoustic potential reflection and transmission coefficients. Tables 9 and 10 show the 
scaled pressure reflection coefficients for nominal flow (left to right) and reverse flow (right to 
left) with the source at x - 0 . These correspond with the left and right hand sides of equation 
(5). The reflection matrices shown in these two tables are seen to be reciprocals, as predicted. 
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mode 

1 

2 

3 

1 

0.16566 + 0.11529 i 

0.23310 + 0.05472 i 

-0.00056 - 0.00099 i 

2 

0.20720 + 0.09263 i 

-0.02814 + 0.01965 i 

0.00035 - 0.00035 i 

3 

-0.00044 - 0.00067 i 

-0.00069 -0.00111 i 

-0.291 (-4) + 0.286(-5) i 


Table 9. Scaled acoustic potential reflection coefficients for direct flow and source at x = 0 for 
the converging circular acoustically treated duct. Corresponds to the left hand side of equation 
(5). r) r = 10 , m = 3. 


mode 

1 

2 

3 

1 

0.16566 + 0.11529 i 

0.20720 + 0.09623 i 

-0.00044 - 0.00067 i 

2 

0.23310 + 0.05472 i 

-0.02814 + 0.01965 i 

-0.00069 - 0.00111 i 

"■> 

J 

-0.00056 - 0.00099 i 

0.00035 - 0.00035 i 

-0.291 (-4) + 0.286(-5)i 


Table 10. Scaled acoustic potential reflection coefficients for reverse flow and source at x = 0 
for the converging circular acoustically treated duct. Corresponds to the right hand side of 
equation (5). t) = 10 , m = 3 . 

Tables 1 1 and 12 show the scaled potential transmission coefficients for nominal flow 
(left to right) with the source plane at x - 0 and reverse flow (right to left) with the source at 
x = L . These correspond with the left and right hand sides of equation (7). The reflection 
matrices shown in these two tables are seen to be reciprocals for modes n = 1 , n = 2 . 
Reciprocity involving mode n - 3 seems to fail. The reason for this can be deduced by 
referring back to Figure 8 and noting that the deeply cut off mode n = 3 creates acoustic 
pressure levels at x = 0 which are probably unresolvable with accuracy by the numerical 
model. To test this hypothesis, another test of reciprocity with the same geometry, flow, and 
mode number was carried out, but with the non-dimensional frequency increased to = 12 
( from r\ = 10 ). This changes the cut-off ratios for the modes n = \,n=2,n=3 to 
5 = 2.88 , % = 1.51 , E, = 1.07 at x = 0 and E, = 1.98 , \ = 1.04 , £ = 0.73 at x = L . 
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This makes three propagating modes at the left end and two propagating modes at the right end, 
and retains the interesting feature of the transition from cut on to cut off for mode n - 3 in a 
transition from left to right, or from cut off to cut on in the opposite direction. 


mode 

1 

2 

3 

1 

-0.18538 + 0.24017 i 

0.04016 + 0.13443 i 

-0.223(-5) - 0.286(-6) i 

2 

0.374(-4) - 0.570(-4) i 

-0.1 15(-4) - 0.288(-4) i 

0.1 10(-8) - 0.344(-8) i 

3 

-0.52 1 (-7) + 0.1 16(-6) i 

0.301 (-7) - 0.507(-7) i 

0.260(-l 1) + 0.595(-l 1) i 


Table 11. Scaled acoustic potential transmission coefficients for direct flow and source at x = 0 
for converging circular acoustically treated duct. Corresponds to the left hand side of equation 
(7). r| r = 10 , m = 3 . 


mode 

1 

2 

3 

1 

-0.18538 + 0.24017 i 

0.374(-4) - 0.570(-4) i 

-0.521(-7) + 0.116(-6)i 

2 

0.04016 + 0.13443 i 

-0.1 15(-4) - 0.288(-4) i 

0.301(-7) + 0.507(-7) i 

3 

0.00014 - 0.00005 i 

0.332(-7) + 0.961(-8) i 

-0.644(-10) - 0.770(-ll) i 


Table 12. Scaled acoustic potential transmission coefficients for reverse flow and source at 
x = L for the converging circular acoustically treated duct. Corresponds to the left hand side of 
equation (7). = 10 , m = 3 . 
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mode 

1 

2 

3 

1 

-0.08862 + 0.41280 i 

0.31074 - 0.05362 i 

0.03674 + 0.06580 i 

2 

-0.02353 - -0.00062 i 

0.02811 - 0.03804 i 

0.00839 - 0.00304 i 

3 

0.877(-5) - 0.791 (-5) i 

0.0255(-5) + 0.237(-4) i 

-0.208(-5) + 0.397(-5) i 


Table 13. Scaled acoustic potential transmission coefficients for direct flow and source at x = 0 
for converging circular acoustically treated duct. Corresponds to the left hand side of equation 
(7). q r = 12 , m = 3 . 


mode 

1 

2 

3 

1 

-0.08862 + 0.41280 i 

-0.02353 - 0.00062 i 

0.877(-5) - 0.791(-5)i 

2 

0.31074 - 0.05362 i 

0.02811 - 0.03804 i 

0.255(-5) +0.237(-4) i 

3 

0.03674 + 0.06580 i 

0.00839 - 0.00304 i 

-0.208(-5) +0.397(-5) i 


Table 14. Scaled acoustic potential transmission coefficients for reverse flow and source at 
x - L for the converging circular acoustically treated duct. Corresponds to the left hand side of 
equation (7). = 12 , m - 3 . 


Figure 9 shows acoustic iso-pressure contours for the case of reverse flow, with the 
source at the right end with n = 3 and r\ r = 12 . The contour levels show that at the left end 
the acoustic pressure levels are substantially higher than those shown in Figure 8, and they are 
more accurately resolved by the modeling scheme. 

Tables 13 and 14 show the scaled potential transmission coefficients for nominal flow 
(left to right) with the source plane at x = 0 and reverse flow (right to left) with the source at 
x = L . The non-dimensional frequency is r\ r = 12. Now reciprocity is satisfied (a reciprocal 
relationship of the scaled transmission matrices) to a high level of accuracy. It is concluded that 
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accurate resolution of the acoustic field of the deeply cut off mode n - 3 . This warns that there 
is a practical limit beyond which reciprocity may not be verifiable for cut off modes. 

A final example considers the converging circular duct profile previously shown in 
Figure 6 with the same steady flow Mach number, but in this case treated as one dimensional. 
Inlet and exit Mach numbers for the nominal flow are M = 0.13 and M = 0.71. Acoustic 
propagation is also taken to be one dimensional at non-dimensional frequency T| r = 1.0 based 
on a reference length which is the radius (of an assumed circular cross section) of the duct at the 
nominal source plane x - 0 . This scaling makes everything consistent with the axially 
symmetric duct formulation. Propagation at the chosen frequency has been modeled by the 
axially symmetric formulation using angular mode m = 0 and radial mode n = 1 . Iso- 
potential contours for the steady flow are shown in Figure 6 and iso-acoustic pressure contours 
are shown in Figure 10 (with longer uniform extensions at the two ends, as compared to Figure 
6). Because of the duct contour and steady flow there is noticeable deviation from true one 
dimensional propagation, however there is no indication of substantial scattering into higher 
order modes, which the one dimensional theory necessarily excludes. 

Table 15 is a s umm ary of the reciprocity and power invariance benchmark tests which are 
available in the one dimensional case. The results of (a) and (b) substantiate the reciprocity 
statement of equation (9) and the results of (c) and (d) substantiate the reciprocity statement of 
equation (11). The power results (e)-(l) substantiate the observations based on power 
considerations that power reflection and transmission coefficients are independent of flow 
direction and source location (equations (13) and (14)). As an indication of the comparison 
between the axially symmetric (3 -Dim) duct model and the one dimensional model, power 
reflection and transmission coefficients for the axially symmetric model are shown also. The 
invariance of the power transmission and reflection coefficients to flow direction and source 
location is true (numerically) at low frequencies in the axially symmetric duct and the one 
dimensional prediction's of reflection and transmission characteristics quite favorably correlate 
with predictions of the axially symmetric model. The properties of invariance of the power 
reflection and transmission coefficients is not generally true at higher frequencies in the axially 
symmetric model when scattering into higher or lower modes occurs. 
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Coefficient 

3-Dim 

1-Dim 

1-Dim 

(a). Reflection coefficient, 
direct flow, source left 



0.430129 + 0.050053 i 

(b). Reflection coefficient, 
reverse flow, source left 



0.430129 + 0.050053 i 

(c). Scaled transmission coefficient, 
direct flow, source left 



0.384607 + 0.145370 i 

(d). Scaled transmission coefficient, 
reverse flow, source left 



0.384607 + 0.145370 i 

(e). Power reflection coefficient, 
direct flow, source left: 

0.171105 

0.187517 


(f). Power reflection coefficient, 
reverse flow, source left 




(g). Power reflection coefficient, 
direct flow, source right 

0.171105 

0.187517 


(h). Power reflection coefficient, 
reverse flow, source right 

0.171105 

0.187517 


(i). Power transmission coefficient, direct 
flow, source left 

0.828895 

0.812483 


(j). Power transmission coefficient, reverse 
flow, source left 

0.828895 

0.812483 


(k). Power transmission coefficient, direct 
flow, source right 

0.828895 

0.812483 


(1). Power transmission coefficient, reverse 
flow, source right 

0.828895 

0.812483 



Table 15. Summary of scattering matrix reciprocity and power invariance benchmark tests for a 
one dimensional converging duct, with power invariance comparisons for an axially symmetric 
model. M q = 0.13, M L = 0.71, r\ r = 1.0 (m = 0 , n = 1 in the axially symmetric case). 
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CONCLUSION 


Numerical verification of the reciprocal relationships derived in References [1,2] has 
been accomplished using a finite element model for duct propagation. Three cases have been 
presented for an axially symmetric duct model, one introducing the feature of transition from an 
annular to a circular duct without and with acoustic treatment, and the second introducing a 
converging duct with substantial steady flow acceleration and segmented acoustic treatment. A 
fourth case is presented for the converging duct at low frequency where a one dimensional 
model of propagation is appropriate. Reciprocal characteristics of the scattering matrices are 
verified with exceptional accuracy, as are predicted relationships for power reflection and 
transmission coefficients in the one dimensional case. 
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Figure 4. Acoustic iso-pressure contours for acoustically treated annular/circular duct with 
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Figure 7. Acoustic iso-pressure contours for acoustically treated converging circular duct with 
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